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HIGH  RESOLUTION  AIR-SEA  INTERACTION  DRI 
GOALS  AND  ACCOMPLISHMENTS 

Technical  goal: 


•  Next  generation  LES  model  of  the  marine  PBL  with  a 
phased  resolved  spectrum  of  surface  waves  A  >  0(1  Ora 

*  Use  measured  2D  wave  fields  as  surface  boundary 
^conditions 

Science  goals: 


•  Wind  response  to  transients,  mis-aligned  winds  an 
waves,  non-equilibrium  winds  and  waves 


TECHNICAL  ACCOMPLISHMENTS 


Built  a  3-D  LES  model  for  a  full  marine  PBL  with  a  spectrum  of  resolved 
time-varying  surface  waves  at  the  lower  boundary 


•  Non-orthogonal,  moving  meshes  that  track  the  time  varying  surface 

•  Co-located  arrangement  of  flow  variables  (compact  stencil) 

•  Can  impose  a  realistic  2-D  wave  spectrum  with  random  phases.  Waves  are 
advanced  in  time  using  the  dispersion  relation 

•  Rough-wall  high  Reynolds  number  boundary  conditions 

•  Accounts  for  the  wave  motion  (orbital  velocities) 

•  Allows  “steep”  waves  ak  ~  0.3 

•  Applicable  to  stationary  terrain  ( e.g 3D  hills  and  bumps) 

•  Mesh  spacing  (Ax,  Ay,  A z)  ~  0(lm)  near  the  water  surface. 

•  Meshes  are  1024  x  1024  x  512  gridpoints 

•  Fully  parallel  MPI  implementation  (tested  with  as  many  as  8192  cores) 
Manuscript: 

A  large  eddy  simulation  model  of  turbulent  atmospheric  boundary  layers  with 
time-dependent  and  stationary  surface- following  grids 


RECENT  LES  CODE  IMPROVEMENTS 


Option  to  impose  Donelan  (1985)  empirical  wave  spectrum 

E(<jj,(j))  =  D(lj,  (j))F((jj) 

Fixes  to  the  surface  boundary  conditions  for  the  wave  orbital  velocities  given 
wave  height  h  =  h(x,y,t) 

dh  dh  dh 

WQ  =  —  +  Uo  - - b  VQ- 


dt 


dx 


dy 


Proper  boundary  condition  for  the  ghost  point  pressure  for  non-orthogonal 
grids  (important  fix  for  waves  and  hills) 

Improved  grid  generation  (smoother  grids) 


Timing  associated  with  advancing  the  wavefields 


•  Ability  to  impose  wave  shapes  from  wind-wave  tank  simulations 

•  Improved  convergence  of  pressure  solver  for  spilling  breakers 
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LES  INVESTIGATIONS  OF  WIND-WAVE  INTERACTION 


LES  database  (incomplete): 

•  Imposed  waves  using  Donelan  spectrum 

•  Neutral  and  weakly  unstable  with  Ug  =[ 5  —  25]  ms-1 

•  Wave  age  Cp/Uio  =  [1.1  —  4.8]  * - N6W 

Differences  between  slow  waves,  fast  waves  and  stationary  bumps: 

•  Flow  structures 

•  Velocity  and  pressure  fields 

•  Statistics  (means,  variances,  fluxes) 

•  Clear  dependence  on  wave  age  in  surface  layer 

•  Non-equilibrium  situations  are  intriguing 

Manuscript: 


Large  eddy  simulations  of  marine  atmospheric  boundary  layers  with  a  resolved 

spectrum  of  surface  waves 


C(m) 


VERTICAL  PROFILE  OF  MEAN  WIND 


WORK  STILL  IN  PROGRESS 


Fine  mesh  resolution  runs  on  “Yellowstone” 
Simulations  of  laboratory  breakers 


Site  simulations  of  HiRES  will  not  be  completed  with  the  current  DRI 
funding 


“YELLOWSTONE"  COMPUTATIONAL  AWARD 


The  NCAR-Wyoming  supercomputer  center  is  home  to  Yellowstone 

•  The  largest  IBM  Idataplex  in  operation 

•  Peta-scale  machine  with  72,000  processors 

Awarded  Proposal:  Large  eddy  simulations  of  high  wind  atmospheric  boundary 
layers  above  a  spectrum  of  resolved  moving  wind  waves 

•  Total  computational  hours:  12  x  10° 

•  Duration:  1  January  2013  -  31  December  2013 

•  Hours  consumed  to  date:  5  x  106 


LES  EXPERIMENTS  ON  YELLOWSTONE 


•  Slightly  unstable  flow  Q *  =  0.01  K  m  s-1 

•  Geostrophic  winds  Ug  =  [5,25]  m  s-1 

•  zi  ~  400  m,  zo  =  0.0002  m 

•  Donelan  wave  spectrum 

•  Wave  age  Cp/Uiq  =  [1.1  —  4.8] 

•  Domain  (3000  x  3000  x  800)  m 

•  Grid  (1024  x  1024  x  512)  « -  0.5  X  109  gridpoints 

•  Ax  =  Ay  =  2.9  m,  Az  =  1.0  m  at  surface 

•  2048  processors,  ~  500,000  “Yellowstone”  core  hours 


Multiple  scales  in  the  MABL:  Can  we  capture  them? 


SURFACE  FLUXES  AND  WIND  SPEED  FROM  HIRES 
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CLOUD  STREETS  OVER  THE  OCEAN 


The  influence  of  wave  breaking  on  the  surface 
pressure  distribution  in  wind-wave  interactions 


By  MICHAEL  L.  BANNER 


FiQi'RR  8.  Photographs  showing  typical  breaking  and  incipient  breaking  waves  in  the  propagating  wind-wave 

lengthscale  is  in  cm.  The  wind  speed  was  5.5 m/s  from  left  to  right,  (a)  3.35  Hz  incipient  breaking  waves.  (5)  3.35  II*  breaking  uaies. 
(r)  2.85  Hz  incipient  breaking  waves.  ( d )  2.85  Hz  breaking  waves. 


FORM  DRAG  OVER  BREAKING  LABORATORY  WAVES 
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SUMMARY 


Coupled  a  LES  for  the  marine  PBL  to  a  realistic  3D  time-dependent 
wavefield 

Coordinate  transform  x*  — >•  &  accounts  for  mesh  movement 

LES  of  neutral  and  weakly  unstable  stratified  flow  above  resolved  waves 

Clear  depdendence  on  wave  age  in  surface  layer 

Competition  between  “wave  pumping”  and  turbulence 

Carrying  forward  high  resolution  simulations  with  109  gridpoints  on 
Yellowstone  using  2048  or  more  cores 

New  LES  for  laboratory  breaking  waves 

HiRES  site  simulations  are  not  planned  under  the  current  DRI 

Two  manuscipts  under  construction 

—  Algorithm  Description 

—  LES  process  studies  of  wind-wave  interaction  for  a  full  marine  PBL 


MEETING  PAPERS  AND  PRESENTATIONS 


19th  Boundary  Layer  &  Turbulence  Meeting,  2010  (2  extended  abstracts) 

20th  Boundary  Layer  &  Turbulence  Meeting,  2012 

Direct  &  Large  Eddy  Simulation  9,  2013 

6th  Theoretical  Fluid  Mechancis  Conference,  AIAA,  2011 

Gulf  of  Mexico  Oil  Spill  and  Ecosystem  Conference,  2013 
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TURBULENT  FLOW  PAST  A  STEEP  3D  HILL 


•  Multi-scale  vortical  wake  flow 


•  Intermittent  pair  of  vortices  near  hill  base 

•  Complex  separation  over  hill  crest  and  flanks 


•  3D  flow  patterns  distinct  from  2D 
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VERTICAL  PROFILE  OF  MEAN  WIND 
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MABL  +  PHASE  RESOLVED  WAVES 


LES  of  a  full  marine  atmospheric  boundary  layer 

-  Rotation,  stratification,  geostrophic  pressure  gradients,  surface  heati 
(cooling) 

3D  time  varying  surface  wave  field 

-  Donelan  et  al(  1985)  spectrum 

-  Measured  spectrum  (WAMOS  vs  ATM)? 

-  Need  for  unresolved  waves  (rough  wall  boundary  condition) 


Funded  computational  proposal  12  million  core  hours  on  “Yellowstone” 

-  10242  x  512  gridpoints,  (Ax,  Ay,Az)  =  (1.17,1.17,0.5)  m 

-  4096  or  8192  processors 


MABL  +  PHASE  RESOLVED  WAVES 


Discussion  points: 

-  Drag  at  high  winds,  viscous  and  pressure  contributions,  separation 

-  Wave  age  regime,  is  wind-wave  equilibrium  the  most  interesting,  non¬ 
equilibrium  regimes  where  waves  may  play  more  of  a  role,  or  growing 
waves 

-  Estimate  of  wind  input  for  spectral  models 

-  Dynamics  of  wind-wave  coupling 

-  Comparison  with  measured  data 


LES  MODEL  EQUATIONS 
FOR  PBL  FLOW  OVER  3D  WAVES  h(x,y,t ) 
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Plus  rough  wall  boundary  conditions  and 
matching  to  orbital  velocity  of  wavefield 


SCALING  OF  FLAT  LES  CODE 


NP  Processors 

Figure  6:  Performance  of  the  NCAR  parallel  LES  code  on  the  NERSC  Cray  supercomputers,  Franklin 
and  Hopper.  NP  is  the  total  number  of  processors  and  the  vertical  axis  is  total  computational  time 
t  x  NP  divided  by  total  work.  Nz  is  the  number  of  vertical  levels  and  Mx>y  is  proportional  to  the  FFT 
work,  i.e.,  MXiV  =  Nx.y\ogNx_y  with  Nx_y  the  number  of  gridpoints  in  the  x  and  y  directions.  Ideal  scaling 
corresponds  to  a  Hat  line  with  increasing  number  of  processors.  This  Hgurc  illustrates  that  the  LES  code 
exhibits  strong  scaling  for  different  combinations  of  problem  size  and  2D  domain  decompositions,  a) 
Green  ♦  problem  size  5123;  b)  Red  ♦  10243;  c)  Black  ♦  20483;  d)  Orange  ♦  30723;  and  c)  Blue  ♦  30723. 
Cases  [a)  -  d)]  arc  timing  tests  performed  on  Franklin  while  ease  c)  is  a  timing  test  on  Hopper.  For  a 
given  number  of  total  processors  NP  the  symbols  arc  different  vertical  and  horizontal  decompositions, 
i.e.,  different  combinations  of  (NPz,NPxy)  as  depicted  in  Fig.  5. 
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ABSTRACT 


e  A  massively  parallel  large-eddy  simulation  (LES)  code  for  atmospheric  planetary  bound- 
7  ary  layers  (PBLs)  above  moving  and  stationary  waves  is  described. 
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1.  Introduction 


9  Large-eddy  simulation  (LES)  is  a  powerful  computational  tool  for  modeling  the  atmo- 

10  spheric  and  oceanic  planetary  boundary  layers  (PBLs)  as  simulations  permit  a  systematic 

11  study  of  stratified  3-D  turbulence  and  its  coupling  with  numerous  small  and  large  scale  pro- 

12  cesses.  Also,  a  key  use  of  simulation  output  is  as  a  virtual  database  to  test,  design  and  build 

13  simpler  1-D  single-column  boundary-layer  parameterizations  for  larger  scale  models  ( e.g ., 

14  Large  et  al.  1994;  Ayotte  et  al.  1995;  Beare  et  ah  2006;  Moeng  et  ah  2010;  McWilliams 

15  et  ah  2012).  Given  this  growing  importance  of  LES  in  studying  boundary- layer  dynamics 
is  it  is  then  critical  to  improve  its  subgrid-scale  parameterizations  as  advocated  by  Wyngaard 
17  (1998,  2004),  but  also  advance  LES  as  a  computational  method  and  further  its  ability  to 
is  couple  atmospheric  (or  oceanic)  turbulence  with  other  processes;  the  work  outlined  here 

19  focuses  on  the  latter. 

20  In  the  past,  LES  of  atmospheric  boundary-layer  flows  was  primarily  restricted  to  ideal- 

21  ized  flows  with  flat  stationary  upper  and  lower  boundaries  owing  to  limitations  in  numerical 

22  schemes  and  computer  power.  Recent  algorithmic  advances  such  as  development  of  efficient 

23  elliptic  solvers,  use  of  non-staggered  variable  layouts,  and  refinements  to  immersed  bound- 

24  ary  techniques  (e.g.,  Zang  et  al.  1994;  Ferziger  and  Peric  2002;  Mittal  and  Iaccarino  2005; 

25  Lundquist  et  al.  2010)  along  with  the  advent  of  massively  parallel  computing  architectures 

26  have  broadened  the  landscape  of  turbulent  flows  within  reach  of  simulations  to  include  at- 

27  mospheric  PBLs  with  boundary  shapes  of  varying  complexity.  Boundary-layer  flows  over 

28  hills  and  moving  water  waves  continue  to  be  usefully  modeled  with  second-order  (or  RANS) 

29  closure  technologies  despite  their  restrictive  assumptions  (e.g.,  Gent  and  Taylor  1976;  Gent 

30  1977;  Chalikov  1978;  Li  1995;  Gong  et  al.  1996;  Taylor  1998;  Belcher  and  Hunt  1998;  Ayotte 

31  and  Hughes  2004;  Chalikov  and  Rainchik  2011).  Past  experience  has  shown  that  the  pre- 

32  dictive  capabilities  of  these  models  are  dependent  on  their  turbulence  closures,  and  in  many 

33  circumstances  the  models  are  linearized  and  thus  only  applicable  to  hills  and  waves  of  small 

34  waveslope  ak;  the  mean  flow  remains  steady  and  fully  attached  under  these  conditions.  The 
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35  shortcomings  of  RANS  closure  models  can  in  principle  be  overcome  by  LES. 

36  The  present  work  follows  similar  paths  as  previous  LES  for  turbulent  flow  over  ide- 

37  alized  rough  stationary  hills  (e.g.,  Gong  et  al.  1996;  Brown  et  al.  2001)  but  is  new  in 

38  that  our  scheme  couples  a  turbulence-resolving  large-eddy  simulation  of  the  atmospheric 

39  boundary  layer  assuming  a  Boussinesq  flow  model  (Sullivan  and  Patton  2011)  to  a  three- 

40  dimensional  time-dependent  resolved  surface  gravity  wavefield,  i.  e. ,  it  is  applicable  to  moving 

41  wavy  boundaries.  The  numerical  scheme  is  based  on  a  curvilinear  coordinate  transformation 

42  from  physical  to  computational  space  with  a  moving  grid.  A  main  motivation  for  constructing 

43  and  applying  this  enhanced  LES  model  is  to  improve  our  understanding  of  air-sea  interaction 

44  dynamics  at  time  and  space  scales  commensurate  with  the  winds,  waves,  and  currents  that 

45  exchange  momentum  and  scalar  fluxes  in  the  lower  atmosphere  and  upper  ocean  (Sullivan 

46  and  McWilliams  2010).  Process  level  understanding  of  the  coupling  between  turbulent  winds 

47  and  waves  in  equilibrium  and  non-equilibrium  situations  is  also  needed  to  further  develop 

48  the  next  generation  of  climate,  weather,  and  wave  forecast  models.  For  example,  coupled 

49  wind- wave-ocean  models  (Chen  et  al.  2007;  Black  et  al.  2007)  are  viewed  as  critical  tools 

50  for  accurate  prediction  of  tropical  cyclone  intensity  and  track  forecasts,  but  these  model- 

51  ing  systems  employ  a  suite  of  parameterizations  that  are  largely  statistical  descriptions  of 

52  the  wind-wave  interactions  that  generate  the  critical  momentum  and  scalar  fluxes.  These 

53  forecast  models  do  not  account  for  the  important  phase  relationships  between  winds,  waves 

54  and  currents,  e.g.,  the  spatial  and  temporal  intermittency  of  wave  breaking  that  occurs  in 

55  moderate  to  high  winds.  Also,  there  is  a  growing  appreciation  that  wave-current  interac- 

56  tions  are  important  for  the  upper  ocean  boundary  layer  at  the  time-scale  of  weather  events 

57  (Sullivan  et  al.  2012),  and  for  climate  predictions  (Belcher  et  al.  2012),  and  that  remotely 

58  generated  swell  and  non-equilibrium  wave  states  can  play  an  important  and  critical  role  in 

59  the  surface-layer  dynamics  of  the  atmospheric  boundary  layer  (Hanley  et  al.  2010;  Sullivan 

60  and  McWilliams  2010). 

61  Our  work  continues  past  efforts  which  coupled  a  turbulent  boundary-layer  flow  with  a 
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62  single  monochromatic  wave  (Sullivan  et  al.  2000;  Sullivan  and  McWilliams  2002;  Sullivan 

63  et  al.  2008;  Nilsson  et  al.  2012),  but  is  an  advance  over  these  earlier  developments.  The 

64  computational  method  described  here  allows  for  nearly  arbitrary  3-D  wavefields,  i.e.,  the  sea 

65  surface  elevation  h  =  h(x,y,t),  as  a  surface  boundary  condition  where  h(x,y,t)  is  assumed 

66  to  be  a  single-valued  function.  The  spatial  scales  of  the  resolved  turbulence  and  waves  are 

67  0( m)  up  to  the  height  of  the  atmospheric  PBL  z%  ~  0(500  m)  or  more.  To  cover  this  broad 
es  range  of  scales  requires  a  large  number  of  mesh  points  and  significant  computational  power. 

69  High  resolution  simulations  of  atmospheric  turbulence  in  the  presence  of  surface  waves  has 

70  the  potential  to  provide  new  insight  into  the  dynamics  of  air-sea  coupling  at  small  scales. 

71  Also,  posing  LES  simulations  in  the  presence  of  both  resolved  and  SGS  waves  (Nakayama 

72  et  al.  2009)  touches  on  similar  parameterization  issues  in  large-scale  atmospheric  models. 

73  In  the  current  problem  posing,  the  surface  waves  are  externally  imposed  based  on  existing 

74  empirical  wave  spectra.  These  spectra  contain  no  information  about  the  phase  of  individual 

75  waves  or  wave  groups  over  a  horizontal  patch  of  the  ocean  surface  that  might  typically  be 

76  covered  by  an  atmospheric  LES  model  (e.g.,  an  area  of  say  3  km  x  3  km).  Theoretical 

77  work  by  Sajjadi  et  al.  (2013)  suggests  that  wave  groupiness  is  potentially  important  for 

78  wave  growth;  this  can  be  partly  accounted  for  in  simulations  by  using  direct  observations 

79  of  the  sea  surface  from  field  campaigns  (Romero  and  Melville  2010).  Although  our  main 
so  computational  target  is  building  a  turbulence  resolving  LES  model  for  studying  wind-wave 
si  interactions,  with  minor  changes  the  algorithm  is  easily  adapted  to  atmospheric  PBLs  with 

82  undulating  (stationary)  landscapes  (e.g.,  Sullivan  et  al.  2010). 

83  The  outline  of  the  paper  is  as  follows:  Section  2  is  an  introduction  to  the  LES  equations 

84  appropriate  for  a  high- Reynolds  number  atmospheric  PBL  and  their  transformation  to  time- 

85  dependent  curvilinear  coordinates;  Section  3  outlines  the  numerical  method;  Sections  4  and 
se  5  describe  typical  applications;  and  Section  6  provides  a  summary  of  the  findings. 
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87 


2.  Governing  equations 


88  In  the  description  of  the  model  equations,  given  below,  the  following  notation  is  used:  u  = 

89  Uj  =  (u,  v,  w)  denote  the  Cartesian  velocity  components,  6  is  virtual  potential  temperature, 

90  p*  is  the  pressure  variable  normalized  by  density  p,  and  e  is  the  subgrid-scale  (SGS)  energy. 

91  Quantities  with  an  overbar  (  )  are  interpreted  as  LES  spatially  filtered  variables. 


92  a.  Atmospheric  LES  model  equations 

93  The  set  of  spatially  filtered  LES  equations  applicable  to  turbulent  flow  in  the  atmospheric 

94  boundary  layer  under  the  Boussinesq  assumption  is  (e.g.,  see  Moeng  1984;  McWilliams  et  al. 

95  1999) 


dujUi  dp 


_  2 eijk<AjUk  +  5a0(O~Oo)  -  ^  (lb) 

at  oxj  oxj  oxj  oxi 

d9_  _  duj9  _  drw 

dt  dxi  dxi 

I  =  “'isf  “  T>iSii  +  PTm  +  w,  f  ■  (ld) 

In  the  above,  O,  =  (0,  0,  fl  sin<h),  with  the  Earth’s  (rotation,  latitude)  denoted  (12,$),  the 
buoyancy  parameter  f3  =  g/0o  where  the  gravitational  acceleration  is  g ,  the  reference  (still 
air)  virtual  potential  temperature  is  0o ,  and  the  large-scale  external  pressure  gradients  nor¬ 
malized  by  density  dV /dx^  drive  the  boundary-layer  winds.  A  Boussinesq  flow  model  requires 
that  p* ,  the  pressure  variable,  be  determined  from  an  elliptic  Poisson  pressure  equation  to 
enforce  mass  conservation;  this  has  significant  impact  on  the  solution  algorithm  (see  Section 
3d).  Ultimately,  the  SGS  momentum  and  temperature  (or  scalar)  fluxes  (r,j, t,Q,  respec¬ 
tively,  require  modeling  in  the  interior  of  the  flow  and  at  the  lower  boundary  (see  Section 
3f)  to  close  the  system  of  equations.  In  the  SGS  turbulence  kinetic  energy  (TKE)  equation 
(Id),  the  time  evolution  of  e  depends  on  right-hand-side  terms  (in  order):  advection  by  the 
resolved  held;  energy  transfer  between  resolved  and  SGS  motions  where  the  resolved  scale 
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in  strain  rate  Stj  =  1/2  ( dui/dxj  +  duj/dxi );  SGS  buoyancy  flux;  Laplacian  diffusion  with  tur- 

112  bulent  eddy  viscosity  vt\  and  viscous  dissipation  S.  Molecular  diffusion  terms  are  neglected 

113  assuming  the  molecular  Reynolds  number  is  high.  A  discussion  of  the  modeling  of  the  SGS 

114  fluxes  is  postponed  to  Section  2d. 

ns  b.  Coordinate  transformation 

ns  We  adapt  our  LES  model  with  a  flat  bottom  to  the  situation  with  a  three-dimensional 
n7  time-dependent  boundary  shape  h  =  h(x,y,t )  by  applying  a  transformation  to  the  physical 
ns  space  coordinates  (x,y,z)  that  maps  them  onto  computational  coordinates  (£,  77,  C)-  The 

119  computational  mesh  in  physical  space  is  surface  following,  non-orthogonal,  and  time  varying. 

120  Vertical  gridlines  are  held  fixed  at  a  particular  (x,  y )  location  on  the  surface  but  the  lines 

121  undergo  vertical  translation  as  a  function  of  time  t,  i.e.,  vertical  gridlines  are  wave  following. 

122  The  transformation  which  obeys  these  constraints  and  maps  the  physical  domain  to  a  flat 

123  computational  domain  x  =>  £  is  the  rule: 

124  T  =  t  ,  f  =f(x)  =  x  >  T)  —  rj(y)  =  y  ,  C  =  C (t,x,y,z)  .  (2) 

125  The  differential  metrics  dxi/d^j  and  dfi/dxj ,  which  are  needed  in  formulating  the  LES  model 

126  in  curvilinear  coordinates,  are  connected  through  the  mapping  transformation.  This  can  be 

127  done  for  general  transformations  (e.g.,  see  Anderson  et  al.  1984,  p.  252-253),  but  here  we 

128  take  advantage  of  our  simplified  terrain  following  grid  (2).  This  leads  to  the  reduced  set  of 

129  non-zero  metric  relationships: 

130  C t  ZfJ  ,  Vy  1  j 

131  J  ,  fy  j  C z  1/ A;  J  (3) 

132  where  J  is  the  Jacobian.  The  time  dependence  of  the  mapping  appears  in  (3)  where  dz/dt  = 

133  zt  is  the  grid  speed,  i.e.,  the  vertical  velocity  of  individual  gridpoints.  Since  our  specific 

134  interest  is  flows  with  wavy  lower  boundaries  we  define  the  rule  in  (2)  that  maps  physical 
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135 


space  onto  flat  computational  space  at  any  time  t  according  to  the  prescription 


136 


z 


C  +  h(x,y,t ) 


(4) 


137  where  is  the  top  of  the  computational  domain  and  w  controls  how  rapidly  the  horizontal 
ns  gridlines  become  level  surfaces  in  physical  space.  Typically  we  use  w  =  3  which  represents 
139  an  acceptable  balance  of  vertical  (squeezing,  stretching)  of  near  surface  grid  cells  at  the  wave 
mo  (crests,  troughs),  respectively.  Note  for  w  —  1,  (4)  reduces  to  £  =  (z  —  h)/(l  —  h/Z^)  which 
mi  is  commonly  used  to  map  stationary  terrain  onto  flat  computational  coordinates  (e.g.,  Henn 
M2  and  Sykes  1999).  (4)  is  a  single  valued  mapping  which  produces  smooth  continuous  metric 
M3  derivatives  (d^i/dxj,dC,i/dt)  in  the  interior  of  the  computational  domain  depending  on  the 
U4  smoothness  of  the  boundary  shape  derivatives  ( ht ,  hXl  hy ). 
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157 


(Thomas  and  Lombard  1979;  Demirdzic  and  Peric  1990)  L  For  our  surface  following  grid 

158  the  GCL  simply  reduces  to  (8b)  since  d^/dt  =  (0,  0,  —zt  J ).  Section  3c  describes  how  (7)  is 

159  used  to  track  the  motion  of  the  grid. 

leo  The  complete  set  of  LES  equations  in  computational  coordinates  under  the  time-dependent 

lei  surface-following  transformation  (2)  and  (3)  is  then: 


162 


163 


164 


165 


166 


167 


m 

dii 

d_  /1\  _  dzt 

dt  d( 

mil)  +  dij  K uj-hjZt)ui] 

§t(j)  +  ^im-***® 

m  0) +  Wj m  ~ Ssi  2,)el 

d  ~l  d &  d^ra  dp*' 
d^i  J  dxj  dxj  d^m 


0 

0 

J 

M 

J 

77 

7 
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(8a) 

(8b) 

(8c) 

(8d) 

(8e) 

(8f) 


lea  (8a)  is  the  mass  conservation  (continuity)  equation,  (8c)  is  the  momentum  transport  equa- 

169  tion,  (8d)  is  the  transport  equation  for  potential  temperature,  (8e)  is  the  subgrid-scale  energy 

170  transport  equation,  and  (8f)  is  the  pressure  Poisson  equation.  The  right  hand  sides  of  (8c, 

171  8d,  8e)  model  physical  processes  in  the  atmospheric  boundary  layer,  viz.,  pressure  gradients, 

172  Coriolis  rotation,  divergence  of  subgrid-scale  fluxes,  buoyancy,  and  in  the  case  of  the  SGS  e 

173  equation  also  turbulent  diffusion  and  viscous  dissipation.  For  completeness  these  terms  are 

1Thomas  and  Lombard  (1979,  see  p.  1032)  develop  the  geometric  conservation  law  using  an  integral 
formulation  but  note  that  it  can  equally  be  derived  from  the  differential  form  of  the  governing  equations  by 
setting  p  =  1  and  iq  =  0,  which  is  the  procedure  used  here. 


174 


gathered  here: 


175 

176 

177 

178 

179 
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-  2e 


ijk 


P*  dfj 
J  dxi 


- 

3  j 


+  £3 i  -Jq- {0  -  0O) 

1  dV 

J  dxi 


d_ 

dti 


d_  f  Tie  d£j  \ 

d V  J  dxi) 

TjjSjj  d  f  2ut  de  d£m  d£j 

J  J  d£j  V  J  dim  dxk  dxk 


Tjk  djj  \ 

J  dxkJ 


£ 

1 


(9a) 

(9b) 

(9c) 


Momentum  and  scalar  advection  are  compactly  written  in  strong  flux- conservation  form 
using  the  volume  flux  or  “contravariant  flux”  velocity  (Rhie  and  Chow  1983;  Zang  et  al. 


181  1994) 

182 


Uj  oi  j 

J  dxj 


(10) 


183  where  Ut  is  normal  to  the  surface  of  constant 

184  The  time  dependence  of  the  grid  modifies  the  LES  equations:  the  Jacobian  appears 

185  inside  the  time  tendency  of  each  transport  equation  and  as  anticipated  vertical  advection 
iso  contains  a  contribution  from  the  vertical  grid  movement,  i.e.,  the  total  vertical  flux  of  variable 
i87  if}  depends  on  the  difference  between  the  contravariant  flux  velocity  and  the  grid  speed 
las  {W  —  zt )7  Examination  of  (8)  shows  that  if  the  velocity  and  scalar  fields  are  set  to  constant 

189  values  then  the  left  hand  sides  of  (8c-e)  reduce  to  (8b).  Hence  the  numerical  method  needs 

190  to  satisfy  the  reduced  form  of  the  GCL  discretely  in  order  to  prevent  artificial  sources  and 

191  sinks  from  developing  in  the  computational  domain,  see  Section  3c. 


192  d.  Parameterization  of  SGS  motions 

193  Formally,  the  LES  equations  are  obtained  by  applying  a  spatial  filter,  term-by-term,  to 

194  the  full  governing  equations  ( e.g Wyngaard  2004).  The  filter  is  assumed  to  be  homogeneous 

195  in  space  and  time,  and  as  a  result  the  operation  order  of  differentiation  and  filtering  can 
we  be  interchanged.  For  boundary  layer  flows  this  assumption  is  violated  and  a  commutation 
197  error  occurs  near  walls  (Ghosal  and  Moin  1995)  -  but  the  error  is  almost  always  ignored 
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in  LES  implementations  as  the  inaccuracies  in  SGS  wall  modeling  are  viewed  as  of  greater 

199  importance  (Sullivan  et  al.  2003).  The  parameterization  problem  is  even  more  pronounced 

200  with  a  multi-mode  moving  lower  boundary  because  of  filtered  products  of  fluctuating  metric 

201  coefficients  and  velocity  (Chalikov  1978),  and  extra  stress-like  terms  arise  from  filtering  the 

202  boundary  shape  (Nakayama  et  al.  2009).  The  parameterization  problem  is  then  similar  in 

203  character  to  representing  SGS  topography  and  canopy  turbulence  in  atmospheric  models; 

204  most  often  their  impact  is  accounted  for  by  a  bulk  drag  parameterization  {e.g.,  Wilson 

205  2002).  In  the  present  application,  like  our  flat-bottom  LES,  the  flow  is  assumed  to  be 

206  laterally  periodic  and  our  wave-following  grid  transformation  results  in  x  —  £  and  y  =  rj. 

207  Thus  we  adopt  the  simplest  approach  and  employ  the  filtering  operation  and  standard  SGS 

208  models  as  in  our  flat  LES  code,  but  account  for  the  varying  Liter  width.  The  LES  equations 

209  are  closed  using  the  SGS  parameterizations  outlined  by  Deardorff  (1972)  and  analyzed  by 

210  Moeng  and  Wyngaard  (1988).  Also,  the  solutions  are  explicitly  filtered,  i.e.,  dealiased,  in 
2n  £  —  r\  planes  at  the  end  of  each  timestep.  The  parameterization  utilizes  an  eddy  viscosity  for 

212  SGS  momentum  and  temperature  fluxes,  and  the  classic  Lilly-model  for  viscous  dissipation: 

213 

d0  e3/2 

214  T"ij  2 Sij  ,  Tjg  ,  £  ( V  ^  —  .  (H) 

215  The  turbulent  eddy  viscosity  ut  =  C*, A^e  and  turbulent  diffusivity  uh  =  (l  +  2£/ A)ut  evolve 

216  under  the  action  of  the  transport  equation  for  SGS  TKE  (8e).  The  constants  (Ck,Cs)  = 

217  (0.1,  0.93)  and  corrections  to  the  stability  corrected  length  scale  £  are  found  in  Moeng  (1984); 

218  Moeng  and  Wyngaard  (1988).  The  Liter  length  scale  A  accounts  for  the  varying  averaging 

219  volume  A3  =  (3/2)2(A£A?/A£/  J)  with  the  “3/2”  factor  accounting  for  dealiasing.  Because 

220  of  movement  and  stretching-squeezing  of  vertical  gridlines  A  varies  with  position  and  time. 

221  To  reduce  the  reliance  on  the  SGS  model  we  use  Lne  resolution  near  the  surface  (Section 

222  4).  Since  the  large  wave  motions  are  fully  resolved  the  total  TKE  at  the  surface  is  often 

223  elevated  and  dominated  by  resolved  variances  further  reducing  the  impact  of  the  SGS  model 

224  compared  to  the  situation  with  Lat  stationary  walls. 
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225  e.  Wavefield  prescription 

226  To  complete  our  LES  model  for  the  marine  atmospheric  PBL  we  need  to  prescribe  the 

227  height  of  the  surface  wavefield  h(x,  y,  t )  in  physical  space  over  a  typical  horizontal  domain  of 

228  the  LES;  this  observational  data  is  generally  unknown.  In  the  absence  of  a  full  description 

229  of  the  wavefield  kinematics  and  dynamics  (e.g.,  Romero  and  Melville  2010),  we  adopt 

230  simplifications  that  allow  us  to  use  existing  information.  We  construct  time  and  space  maps 

231  of  h(x,  y,  t )  based  on  typical  empirical  fits  of  measured  two-dimensional  wave  spectra 

232  E(k,<f>)  =  S(k)  D(k ,0)  (12) 

233  where  the  amplitude  S(k)  and  directional  D(k,4> )  spectra  depend  on  the  magnitude  of  the 

234  horizontal  wavenumber  k  =  |k|  =  \kxi  +  ky j|  and  wave  direction  <fi.  Note  these  spectra 

235  are  statistical  averages  that  contain  no  phase  information  about  the  wavefield.  For  the 

236  amplitude  and  directional  spectra  we  choose  the  well  accepted  functional  forms  proposed 

237  by  Donelan  et  al.  (1985)  (see  also  p.  187  of  Komen  et  al.  (1994))  which  are  applicable  to 

238  both  equilibrium  and  non-equilibrium  wind-wave  states.  These  empirical  fits  depend  on  bulk 

239  environmental  parameters,  viz.,  the  reference  surface  wind  speed  at  10  m,  U\q,  the  phase 

240  speed  of  the  peak  in  the  wave  height  spectrum  Cp,  their  ratio,  the  wave  age  Cp/Uw,  and  the 

241  mean  wave  propagation  direction  (0).  To  use  these  measured  spectra,  which  are  expressed 

242  in  terms  of  radial  frequency  u>,  in  the  LES  they  are  transformed  into  wavenumber  space 

243  using  the  linear  dispersion  relation  k  =  u2  /  g  and  properly  weighted  so  that  the  resulting 

244  wavenumber  spectra  are  variance  preserving  (Phillips  1977,  p.  105).  In  physical  space,  the 

245  rule  for  constructing  the  synthetic  wavefield  h(x,  y,  t )  is  as  a  sum  of  linear  plane  waves  each 

246  with  random  phase  p 

247  /i(x,  t)  =  A (k)  exp  [i(k  •  x  —  u(k)  t  +  <p)\  ;  (13) 

k 

248  the  amplitude  A2{k)  =  S(k)D(k,(j))dkxdky.  (13)  is  implemented  in  the  LES  code.  Given 

249  an  initial  state  h(x,y,t0),  future  values  of  h(x,y,t )  are  then  efficiently  obtained  by  using 
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250  2-D  Fast  Fourier  Transforms  (FFTs).  Time  and  space  derivatives  of  the  boundary  shape 

251  (. ht ,  hx,  hy ),  which  are  needed  to  construct  the  moving  mesh  and  the  surface  boundary  con- 

252  ditions,  are  also  computed  in  straightforward  fashion  from  (13)  using  FFTs.  We  anticipate 

253  future  LES  applications  will  directly  ingest  measured  maps  of  h(x,  y,  t )  from  airborne  irn- 

254  agery  (e.g.,  Romero  and  Melville  2010)  or  ship  mounted  X-band  radars. 

255  f.  Velocity  boundary  condition 

256  We  need  to  specify  boundary  conditions  on  the  winds  that  are  consistent  with  the  imposed 

257  surface  wavefield.  First,  both  the  wind  and  wave  fields  are  assumed  to  be  spatially  periodic 

258  in  computational  £  —  rj  planes  and  thus  the  important  boundary  conditions  that  need  to  be 

259  set  are  at  the  bottom  and  top  of  the  computational  domain.  At  the  wave  surface,  the  total 

260  time  rate  of  change  of  the  wave  height  is 

dk  r  1  r  ,\ 

261  —  =  Wo  =  ht  +u0hx  +  vQhy  (14) 

262  where  (u0,  vQ,  wQ )  are  the  water  motions.  Initially  we  assume  these  water  motions  are  domi- 

263  nated  by  the  wave  orbital  velocities  {e.g.,  Lighthill  1978),  i.e.,  there  is  no  net  surface  drift. 

264  The  importance  of  surface  drift  for  wave  breaking  is  discussed  by  Banner  and  Melville  (1976); 

265  Banner  (1990).  The  contravariant  velocity  component  normal  to  the  water  surface  is  given 

266  by  (10): 

267  W  =  Uy  +  v*f  +  w  .  (15) 

268  Matching  (14)  and  (15)  and  requiring  no  flow  normal  (across)  the  wavy  boundary  leads  to 

269  W  —  ht  =  0  at  C  =  0  .  (16) 

270  (16)  effectively  implies  that  spray  and  bubbles  do  not  pass  through  the  interface.  At  the 

271  upper  boundary  of  the  computational  domain  the  horizontal  gridlines  are  level  surfaces  in 

272  physical  space  with  =  zv  =  0.  Then 

273  W  =  w  =  zt  at  C  —  %l  1  (17) 
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274  so  that  the  resolved  vertical  flux  of  momentum  and  scalars  is  zero  across  the  top  boundary. 

275  The  specification  of  the  surface  fluxes  appropriate  for  a  high  Reynolds  number  rough  wall 

276  LES  model  taking  into  account  the  wave  motion  is  described  in  Section  3f  with  additional 

277  details  provided  in  the  Appendix. 

278  3.  Numerical  method 

279  a.  Spatial  discretization  and  variable  layout 

280  The  LES  set  of  equations  (8)  is  discretized  using  the  well-developed  techniques  in  our  flat 

281  bottom  boundary  codes  (Moeng  1984;  Sullivan  et  al.  1996).  Spatial  derivatives  in  the  (£,77) 

282  computational  coordinates  are  estimated  using  pseudospectral  approximations  while  centered 

283  second-order  accurate  finite  difference  approximations  are  employed  for  (  derivatives.  In 

284  order  to  be  compliant  with  (8b),  the  spatial  differencing  evaluates  the  advective  terms  in 

285  (8c-e)  in  flux  form  as  opposed  to  the  rotational  form  used  by  Moeng  (1984)  or  skew-symmetric 

286  form  used  by  Sullivan  et  al.  (2000). 

287  The  variable  layout  in  the  computational  mesh  uses  a  co-located  arrangement  for  the 

288  fundamental  solution  variables  ( Ui,p*,9,e )  which  is  advantageous  as  it  results  in  a  compact 

289  differencing  stencil  and  is  also  applicable  to  meshes  with  large  bends  in  the  coordinate  lines 

290  (Demirdzic  and  Peric  1990;  Zang  et  al.  1994).  The  fundamental  difficulty  of  tightly  coupling 

291  the  velocity  and  pressure  in  a  co-located  mesh  is  circumvented  by  locating  the  contravariant 

292  flux  velocities  at  cell  faces  mimicking  the  usual  arrangement  in  a  staggered  mesh.  To  satisfy 

293  (8a),  written  for  contravariant  velocities,  with  our  mixed  pseudospectral  finite-difference 

294  scheme  results  in  a  novel  layout  of  flow  variables  with  (U,  V)  also  located  at  cell  centers  and 

295  W  naturally  located  at  the  upper  and  lower  cell  faces  as  shown  in  Fig.  1.  Furthermore, 

296  with  a  surface-following  non-orthogonal  grid  (U,V)  =  ( u,v)/J  are  simply  scaled  versions 

297  of  their  Cartesian  counterparts,  see  (10).  Other  mesh  types  and  differencing  schemes,  e.g., 

298  orthogonal  meshes  and  finite  difference  methods,  result  in  different  couplings  between  the 
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299  Cartesian  and  contravariant  velocities  (e.g.,  Zang  et  al.  1994;  Sullivan  et  al.  2000).  Lastly, 

300  the  flow  is  assumed  to  be  spatially  periodic  in  horizontal  £  —  rj  planes  so  as  to  be  consistent 

301  with  a  Fourier  representation. 


302  b.  Time  advancement 


303  Time  integration  is  a  fully  explicit  third-order  Runge-Kutta  (RK3)  scheme  (Spalart  et  al. 

304  1991;  Sullivan  et  al.  1996,  2008)  that  uses  dynamic  time  stepping  with  a  fixed  Courant- 

305  Fredrichs-Lewy  (CFL)  number  and  employs  a  fractional  step  method  to  enforce  the  diver- 

306  gence  free  condition.  As  the  scheme  is  completely  explicit  it  is  equivalent  to  the  so-called 

307  “pressure  correction”  method  of  Armfield  and  Street  (2002).  The  general  rule  for  advancing 

308  a  cell-centered  Cartesian  velocity  variable  to  a  new  time  level  n  over  a  timestep  At  is 


309 


Ui 


n—  1 


j 


Atoin 


d_ 

d£j 


n 


5 


(18) 


310  where  the  pressure  gradient  is  written  in  conservative  form.  The  intermediate  velocity 


Ui 

J 


n—  1 


Hi 

7 


n—  1 


+  A  tan  — 


n—  1 


+  A 0n  ^ 


n— 2 


(19) 


312  is  the  discrete  sum  of  the  full  right  hand  side  Qi/ J  from  the  previous  time  (or  stage)  step 

313  (n  —  2)  and  the  partial  right  hand  side  qi/  J  from  the  current  time  (n  —  1)  minus  the  pressure 

314  contributions;  (an,  j3n)  are  weights  associated  with  the  RK3  scheme  (see  Sullivan  et  al. 

315  2000).  The  time  integration  rule  for  scalar  fields  is  identical  to  (18)  and  (19)  in  the  absence 

316  of  pressure  gradients.  In  our  LES  equation  set,  the  divergence  free  condition  is  satisfied  at 

317  every  RK3  stage  as  opposed  to  the  method  advocated  by  Le  and  Moin  (1991)  (see  also  Zheng 

318  and  Petzold  2006)  which  projects  the  velocity  onto  a  divergence  free  field  only  at  the  final 

319  stage.  The  latter  scheme  is  computationally  more  efficient  since  it  reduces  pressure  iterations 

320  over  a  full  timestep,  but  both  flux  and  skew-symmetric  forms  of  the  advection  term  in  our 

321  scalar  equations  are  based  on  on  satisfying  (8a)  at  all  stages  of  the  RK3  time  advancement. 

322  Attempts  to  use  the  method  of  Le  and  Moin  (1991)  resulted  in  numerical  instabilities. 
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323 


c.  Geometric  conservation  law 


324  The  geometric  conservation  law  (8b)  needs  to  be  satisfied  using  the  same  time  advance- 

325  ment  scheme  as  for  velocity  and  scalar  fields  (see  (18)).  Given  an  externally  imposed  wave- 

326  field,  h(x,  y,  t )  is  known  at  current  and  future  time  steps  (n  —  1,  n)  and  we  can  then  construct 

327  the  distribution  of  vertical  grid  points  at  “W— points”,  viz.,  z(x,y) |(n_1,n)  according  to  the 

328  map  (4).  Jacobians  J|(n-1>n)  are  then  naturally  built  from  from  the  transformation  rule 

329  (3).  Since  we  choose  to  specify  the  grid  point  locations  at  future  timesteps,  we  must  then 

330  find  matching  grid  speeds  zt(x,y )  so  that  (8b)  is  obeyed  discretely.  The  companion  to  (18), 

331  appropriately  inverted,  for  determining  grid  speeds  is  the  semi-discrete  relationship 


332 


dzt 

d( 


1 

an  At 


-  A  tpn 


(20) 


333  This  first-order  equation  is  integrated  vertically  from  the  surface  to  the  top  of  the  domain 

334  with  the  boundary  condition  zt  —  ht  at  (  —  0.  We  emphasize  that  the  upper  boundary  of 

335  our  computational  domain  is  far  from  the  lower  surface  and  the  wave-following  gridlines  in 

336  the  interior  of  the  domain  quickly  asymptote  to  flat  level  surfaces  with  increasing  (.  As  a 

337  result,  zt  at  a  particular  x  —  y  gridpoint  must  vary  with  distance  from  the  wave  surface. 

338  This  is  different  than  the  scheme  proposed  by  Chalikov  (1998)  where  the  grid  speeds  at  all  ( 

339  are  equal  to  the  vertical  motion  of  the  wavefield,  i.e.,  zt  is  constant  with  height.  The  latter 

340  scheme  also  implies  that  the  computational  mesh  at  any  Q,  including  the  upper  boundary, 

341  mimics  the  shape  of  the  lower  surface  which  complicates  posing  upper  boundary  conditions. 

342  There  is  a  subtlety  in  the  implementation  of  (20).  Since  our  time  stepping  is  explicit  the 

343  value  of  Zt  diagnosed  from  (20)  is  at  time  level  n  —  1  but  uses  information  at  time  level  n. 

344  In  other  words,  if  the  wavefield  is  imposed,  to  diagnose  the  grid  speed  at  time  n  requires 

345  knowledge  about  the  wave  shape  at  a  future  time  n  +  1 .  The  end  result  is  that  the  grid  must 

346  be  stored  at  three  time  levels  in  the  LES  code.  Information  about  the  wavefield  at  each  RK3 

347  stage  is  used  in  two-ways  in  the  grid  generation:  h(x,  y,  t )  is  used  to  position  the  grid  nodes 

348  in  the  mesh  and  ht(x,  y,  t)  is  used  to  determine  how  rapidly  the  grid  nodes  in  the  mesh  move 
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349 


to  their  new  locations. 


350  d.  Pressure  Poisson  equation 


The  crux  of  the  numerical  scheme  is  the  formulation  and  solution  of  a  pressure  Poisson 
equation  that  results  in  a  flow  held  consistent  with  our  incompressible  Boussinesq  how  model. 
Compared  to  a  hat  wall  LES  code,  the  solution  of  the  elliptic  pressure  equation  with  a  wavy 
boundary  is  computationally  more  expensive,  a  consequence  of  the  variable  coefficients  and 
mixed  spatial  derivatives  in  (8f).  As  is  standard  practice  we  develop  a  pressure  equation  by 
combining  the  continuity  equation  and  our  time  stepping  scheme.  Substitution  of  (18)  into 
(10)  leads  to  the  time  advance  for  the  contravariant  velocity 


ur  =  Ui 


A  tlfrn  8 &  d irn  dp* 

J  dxj  dxj  d£m 


where  now  the  pressure  gradient  is  expanded  in  non-conservative  form.  (21)  is  designed  to 
mimic  the  usual  staggered  arrangement  of  how  variables  (in  our  hat  LES  code)  with  the 
pressure  located  at  a  cell  center  for  all  three  velocity  components.  The  key  step  in  imple¬ 
menting  (21)  is  the  construction  of  the  intermediate  velocity  held  Ut.  We  use  momentum 
interpolation  of  the  Cartesian  velocity  components  to  build  the  right  hand  side  of  the  con¬ 
travariant  hux  velocities  as  proposed  by  Zang  et  al.  (1994)  and  Sullivan  et  al.  (2000,  2008). 
Essentially,  we  interpolate  u,JJ  to  the  same  (staggered)  locations  as  Ui  and  thus 

6i  =  §[  Sy  <22> 

where  (  )|j  denotes  an  interpolated  value.  As  mentioned  previously,  because  of  the  spatial 
differencing  scheme  only  W |n_1,  located  at  the  upper  and  lower  cell  faces  is  obtained  by 
interpolation,  see  Section  3a.  The  divergence  of  (21)  leads  to  the  pressure  Poisson  equation 


d  A  d^i_  dtUn  dp* 
d^i  J  dxj  dxj  <9£m 


1  dUj 
d^i 


Our  strategy  for  folding  the  pressure  utilizes  an  iterative  method  that  avoids  directly 
forming  the  complicated  left  hand  side  of  (23).  Instead,  examination  of  (18),  (21),  and  (23) 
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suggests  the  equivalent  stationary  iteration  scheme  to  find  p*\n 

1  dUk(p* 


£(p* ,+1)  =  £{p* 1 )  — 


(24) 


A  tan  d^k 

where  the  source  term  is  the  divergence  free  condition,  superscript  i  denotes  the  iteration 
index,  and  the  linear  preconditioning  operator 


1  (d  V 

£(P*)  =  77V 


+ 


d  2p* 


d  (  .  .  dp* 

+  ac  sc 


(25) 


(J)  \  d£2  '  drf 

Note  Uk,  not  Uk,  appears  in  the  source  term  of  (24).  At  each  iteration  we  update  Uk  using 
the  latest  p*  ‘  field  according  to  (21)  and  then  compute  the  new  right  hand  side  of  (24).  C  is 
designed  to  be  an  easily  invertible  diagonal  approximation  of  the  left  hand  side  of  (23)  with 
the  spatially  averaged  Jacobian 


382  (J)( C)  =l[  C)dfd77.  (26) 

U  A 

383  Since  (J)  is  solely  a  function  of  £  this  allows  (24)  to  be  solved  for  p*  ,+1  using  standard  rneth- 

384  ods,  viz.,  2-D  Fourier  de-composition  in  £  —  p  planes  followed  by  tridiagonal  matrix  inversion 

385  in  the  £  direction.  At  convergence  C(p*  1+1 )  =  C(p*  1 )  and  the  divergence  free  condition  (8a) 

386  is  then  satisfied  at  the  new  time  step  n.  Inversion  of  (24)  for  p* 1+1  is  straightforward  and  can 

387  be  done  in  place,  but  is  expensive  since  its  computational  steps  require  global  communication 

388  in  a  parallel  algorithm  (Sullivan  and  Patton  2011).  Solving  for  pressure  in  the  presence  of 

389  wavy  boundaries  increases  the  overall  computational  cost  of  simulations  by  ~  50%  compared 

390  to  simulations  with  flat  lower  boundaries  where  the  pressure  is  solved  for  directly. 


391  e.  Pressure  boundary  conditions  at  a  wavy  boundary 

392  Compatible  velocity  and  pressure  boundary  conditions  are  naturally  built  into  the  it- 

393  eration  scheme  given  by  (24).  First,  the  appropriate  surface  boundary  condition  on  the 

394  diagonal  preconditioner  C  is  simply  dp*/d(  =  0,  which  is  the  condition  most  often  used  in 

395  flat  wall  LES.  The  pressure  boundary  condition  used  in  the  source  term  dUk(p*) / d£k  of  (24) 
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is  exposed  by  examining  the  update  for  (U,  V,  W)  in  advancing  from  n  —  1  to  n  given  by 

(21): 


u  p  =  ur1 
v\n  =  f/|n_1 


A  tar 
A  tar 


1  dp*  (x  dp 
J  d^  + 

1  dp* 


+ 


J  dC_ 

CydV 


[J  dp  J  d( 


W \n  =  W I”-1  -  A tar 


rCx  +  Cy  +  Cz  dp 


2  Q-T 


Cxdp*  C ydp* 

I  r\  I 


(27a) 

(27b) 

(27c) 


J  d(  J  <9£  J  dr) 

In  the  interior  of  the  computational  domain  the  pressure  gradient  dp*/dC,  which  appears 
in  all  three  equations,  is  evaluated  using  standard  second-order  centered  finite-difference 
formulas,  while  the  pressure  gradients  (dp* / dC ,  dp* / dr))  are  evaluated  using  pseudopsectral 
differencing.  Inspection  of  (27a, b)  shows  that  to  update  (U,  V)  using  a  centered  formula 
for  dp*/d(  at  the  first  level  (  =  A£/2  above  the  wavy  boundary  requires  the  pressure  field 
at  C  =  3A<C/2  and  at  the  ghost  point  (  =  —A(/2  below  the  wavy  boundary,  assuming 
equally  spaced  points  in  (,  see  Fig.  2.  Furthermore  at  (  =  —  A£/2,  U  and  V,  are  globally 
connected  to  an  entire  plane  of  ghost  point  pressures  p*(x,y ,  —  A£/2)  through  the  gradients 
(dp* /<9£,  dp* /dp). 

In  order  to  find  the  proper  ghost  point  pressures  we  use  the  W  equation,  its  surface 
boundary  condition,  and  our  iteration  scheme.  At  C,  —  0,  (27c)  is  rearranged  to  yield 


C  +  Cy  +  Cz  P*(C,  V,  AC/2)  -  P*(C,  V,  -  AC/2) 


J 


AC 


W-ht 
A  tan 


C,xdp*  C ydp* 


[  J  J  drj  _ 


i— 1 


(28) 

where  we  impose  the  boundary  condition  (16).  The  superscript  i  in(28)  is  again  the  iteration 
index,  and  for  clarity  we  drop  the  superscripts  n  and  n  —  1  denoting  the  time  step.  At  each 
iteration  with  a  known  held  p*(C,  r),  AC/2)*,  (28)  is  rearranged  to  solve  for  the  ghost  point 
pressure  p*(C,  r),  —  AC/2)*  with  the  right  hand  side  (centered  on  C  =  0)  lagged  from  the 
previous  iteration.  Subsequently,  (27a, b)  are  updated  at  all  interior  nodes  and  a  pressure 
iteration  is  completed  by  sweeping  through  (24)  with  a  new  right  hand  side.  W( C  =  0), 
which  appears  in  (28),  is  constant  during  pressure  iterations  and  is  obtained  from  interior 
node  interpolations. 
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421  Preliminary  tests  indicated  that  extrapolation  of  the  interior  pressures  to  £  =  0  or  use 

422  of  d p*/d(  =  0  at  the  wave  surface  resulted  in  spoiled  solutions  with  oscillations  in  the 

423  pressure  and  vertical  velocity  fields  propagating  from  the  wavy  surface  into  the  interior  of 

424  the  computational  domain.  The  magnitude  of  the  oscillations  increased  with  finer  resolution 

425  and  with  the  number  of  wave  modes  in  the  boundary  shape.  The  use  of  (28)  permits  us 

426  to  impose  the  maximum  number  of  wave  modes  in  h(x,  y,  t )  allowed  by  the  horizontal  grid 

427  resolution  in  our  pseudospectral  differencing  scheme,  i.e.,  2/3(NX,  NY)  where  (NX,  NY') 

428  are  the  number  of  gridpoints  in  the  (£,  rf)  directions,  respectively.  For  simulations  with  flat 

429  boundaries  the  updates  for  ( U ,  V)  are  independent  of  the  ghost  point  pressures,  while  in  a 

430  fully  coupled  air-water  simulation  the  surface  pressure  would  be  consistently  determined  by 

431  matching  with  the  water  motions. 

432  /.  Surface  fluxes 

433  In  high  winds,  the  ocean  surface  is  assumed  to  be  in  a  fully  rough  regime  with  negligible 

434  contributions  from  the  thin  molecular  viscous  sublayer  (Donclan  1998).  This  assumption 

435  however  needs  to  be  re-considered  at  low  winds  with  very  small  scale  waves  (e.g.,  Harris 

436  et  al.  1996;  Meirink  and  Makin  2001)  where  the  roughness  Reynolds  number  can  be  in  a 

437  smooth  or  transitional  regime.  Here,  we  follow  the  approach  of  second-order  closure  modeling 

438  for  turbulent  flow  over  hills  and  waves  and  assume  the  water  surface  is  fully  rough  so  that 

439  the  surface  fluxes  of  momentum  and  scalars  can  be  represented  in  terms  of  standard  law- 

440  of-the-wall  relationships  (e.g.,  Gent  and  Taylor  1976;  Makin  et  al.  1995;  Belcher  and  Hunt 

441  1998).  This  approach  is  less  certain  in  LES  modeling  where  surface  undulations  (or  waves) 

442  are  partially  resolved,  i.e.,  the  velocity  and  pressure  fields  induced  by  the  large  energetic 

443  components  of  the  waveheld  are  resolved  and  the  smaller  scale  less  energetic  waves  are 

444  unresolved  (i.e.,  subgrid-scale).  How  to  parameterize  surface  fluxes  in  a  mixed  resolved-SGS 

445  regime  for  stationary  waves  is  a  current  research  topic  (e.g.,  Nakayama  et  al.  2004;  Anderson 

446  and  Meneveau  2010),  and  is  unexplored  for  moving  waves.  To  make  progress,  we  adopt  a 
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447  simplified  approach  that  applies  the  bulk  aerodynamic  formulas  point-by-point  along  the 

448  wave  surface  in  a  local  Cartesian  coordinate  system: 

449  tq  =  —ul  =  —  Cd\us\2  (29a) 

450  Q*  {9 s  Qsurf')  •  (29b) 

451  where  (Cd,  Ch )  are  bulk  transfer  coefficients  for  momentum  and  temperature,  w*  is  the  friction 

452  velocity,  9S  is  the  potential  temperature  at  the  first  grid  point  (  =  A(/2  off  the  surface,  9surj 

453  is  the  temperature  of  the  water  surface,  and  Q*  is  the  surface  temperature  (heat)  flux.  us  is 

454  the  vector  difference  between  the  surface  following  winds  at  Q  =  A(/2  and  the  water  motions 

455  at  C  —  0-  The  total  stress  r0  in  surface  following  coordinates  is  converted  into  a  surface  flux 

456  matrix  in  the  LES  Cartesian  frame  of  reference  using  the  transformation: 

457  Tjj  =  T  Ta  bij  ,  (30) 

458  where  the  transformation  matrix  bij  varies  with  each  surface  grid  point.  The  lengthy  de- 

459  tails  of  defining  the  local  wave-fitted  Cartesian  coordinate  system,  computing  (Cd,Ch),  arid 

460  generating  btj  are  confined  to  the  Appendix. 

461  4.  Simulations  with  wavy  lower  boundaries 

462  To  demonstrate  the  capabilites  of  our  LES  method,  a  series  of  simulations  with  varying 

463  geostrophic  wind  Ug  =  (5,10,20,25)  m  s^1  are  carried  out  for  a  weakly  unstable  marine 

464  PBL  in  a  domain  (XL,  YL,  ZL)  =  (3000,3000,800)  m  using  (Nx,  Ny,  Nz)  =  (1024,1024,512) 

465  gridpoints;  thus  the  horizontal  grid  spacing  Ax  =  Ay  =  2.93  m  with  the  first  vertical  level 

466  nominally  located  1  m  above  the  water  surface.  A  slice  of  the  computational  mesh  is  given 

467  in  figure  5.  The  initial  temperature  sounding  9  =  300  K  up  to  the  inversion  height  zt  =  400 

468  m,  beyond  this  height  9  increases  linearly  at  3  x  10-3  K  m_1.  The  surface  heating  Q*  =  0.03 

469  Km  s^1,  the  surface  roughness  za  =  0.0002  m,  and  the  Coriolis  parameter  /  =  10-4  s^1. 

470  The  wavefield  is  built,  as  discussed  in  Section  e,  based  on  a  wind  speed  of  15  m  s'1  and  the 


20 


471  phase  speed  of  the  peak  in  the  spectrum  Cp  ~  18  m  s”1.  Thus  the  suite  of  simulations  allows 

472  us  to  examine  a  wide  variation  of  wave  age  from  swell  dominated  to  wind-wave  equilibrium. 

473  Table  1  lists  bulk  properties  of  the  simulations,  viz.,  the  geostrophic  wind,  wave  age,  and 

474  friction  velocity  w*.  U\o  is  the  reference  wind  speed  at  a  height  of  10  m.  The  simulations  are 

475  run  for  more  than  50,000  timesteps  using  restart  volumes  with  fully  developed  turbulence. 

476  The  iteration  count  in  the  pressure  Poisson  solver  set  to  30  and  the  calculations  run  on  2048 

477  computational  cores. 

478  The  wave  height  held  varies  with  space  and  time  as  indicated  by  (13),  see  Fig.  3.  Closer 

479  inspection  of  the  image  shows  a  preference  towards  long  crested  waves,  i.e.,  waves  with 

480  aspect  ratio  long  in  the  direction  perpendicular  to  the  winds.  Horizontal  spectra  of  the  wave 

481  height  variance  in  the  x  and  y  directions  are  given  in  Fig.  4.  These  are  obtained  from  the 

482  Donelan  et  al.  (1985)  empirical  fits  and  converted  into  wavenumber  space  using  the  linear 

483  dispersion  relation.  The  spectra  show  that  there  is  a  preference  for  generating  long  crested 

484  waves  in  the  direction  perpendicular  to  the  mean  wind  direction  at  low  wavenumber.  The 

485  spectra  exhbit  a  power  law  of  k~xxx  at  small  scales  (high  wavenumbers). 

486  5.  Results 

487  Previous  held  observations  (Grachev  and  Fairall  2001;  Smedman  et  al.  1999),  turbulence 

488  closure  modeling  (Hanley  and  Belcher  2008;  Makin  2008),  and  our  own  idealized  LES  (Sulli- 

489  van  et  al.  2008)  all  show  that  fast  moving  swell  can  induce  marked  changes  in  the  atmospheric 

490  surface  layer  winds,  viz.,  an  upward  momentum  hux  from  the  ocean  to  the  atmosphere,  a 

491  low-level  wind  maximum,  and  departures  from  law-of-the-wall  scaling.  The  preliminary  LES 

492  computations  performed  here  over  a  more  realistic  sea  surface  are  in  good  qualitative  agree- 

493  rnent  with  the  previous  studies  but  suggest  the  impact  of  swell  on  the  surface  layer  winds  is 

494  sensitive  to  the  content  of  the  wave  spectrum. 

495  One  of  the  surprising  results  from  the  present  simulations  is  the  significant  impact  of  swell 
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496  on  the  coherence  and  magnitude  of  the  near-surface  pressure  fluctuations.  This  is  illustrated 

497  in  figure  6  where  we  compare  p'/p  for  two  levels  of  wind  forcing  Ug  =  (5,20)  m  s-1,  i.e., 

498  a  low-wind  situation  with  swell  and  a  high  wind  case  approaching  wind-wave  equilibrium. 

499  The  difference  in  the  pressure  signals  is  striking  and  even  more  remarkable  in  animations 

500  of  the  pressure  held.  In  the  low-wind  swell  case  there  is  a  very  strong  correlation  between 

501  p'/p  <  0  and  wave  crests  and  similarly  between  p'/p  >  0  and  wave  troughs  that  extends 

502  over  the  depth  of  the  surface  layer.  Inspection  of  the  how  visualization  and  animations 

503  reveals  that  the  strong  correlation  persists  across  the  range  of  resolved  waves,  i.e.,  both 

504  large  and  small  scale  waves  appear  to  induce  a  similar  pressure  pattern.  The  coherence  of 

505  the  wave  induced  pressure  held  can  extend  to  20  m  or  more  depending  on  the  amplitude  of 

506  the  underlying  wave.  Also,  the  pressure  signatures  propagate  at  the  speed  of  the  waveheld, 

507  additional  evidence  that  the  signals  are  generated  by  surface  waves  and  not  atmospheric 

508  processes.  These  are  clear  signatures  of  “wave  pumping”  by  the  surface  waveheld  on  the 

509  atmosphere.  The  amplitude  of  the  wave  spectrum  (and  hence  the  level  of  wave  forcing)  is 

510  held  constant  in  our  simulations  but  the  magnitude  of  the  turbulence,  as  measured  by  w*, 

511  increases  substantially  with  increasing  wind  speed.  The  structure  of  the  near  surface  pressure 

512  held  is  a  result  of  these  two  competing  effects.  At  low  winds  coherent  pressure  signals  are 

513  generated  by  the  wave  motions  when  the  turbulence  is  weak  but  this  coherence  is  destroyed 

514  by  strong  turbulence  at  higher  winds. 

515  Figure  7  shows  that  the  impact  of  wave  age  also  appears  in  the  vertical  velocity  helds. 

516  In  the  low-wind  swell  regime  we  observe  large- amplitude  large-scale  fluctuations  in  w'.  At 

517  higher  winds  the  spatial  coherence  of  w'  is  destroyed  by  strong  turbulence.  Note  each  panel 

518  in  figure  7  is  sampled  at  the  sample  height  above  the  waveheld.  Also,  the  helds  are  made 

519  dimensionless  by  friction  velocity  w*  which  further  illustrates  the  strong  impact  of  the  wave 

520  motions  on  the  winds  in  the  surface  layer. 

521  In  figure  8  we  compare  vertical  prohles  of  the  mean  wind  speed  and  turbulence  variances 

522  for  the  different  simulations.  These  statistics  are  computed  by  averaging  in  computational 
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523  coordinates,  i.e.,  across  horizontal  planes  at  constant  vertical  height  £.  Similar  to  Sullivan 

524  et  al.  (2008)  we  find  that  the  wind  speed  and  turbulence  variances  depend  on  wave  age.  At 

525  high  winds  as  the  simulations  approach  wind-wave  equilibrium,  the  non-dimensional  wind 

526  profile  ( U)/u *  smoothly  approaches  the  variation  predicted  by  law-of-the-wall.  Significant 

527  differences  are  observed  for  the  cases  dominated  by  swell:  the  surface  layer  winds  are  acceler- 

528  ated  compared  to  rough  wall  scaling.  As  suggested  by  the  flow  visualization,  the  turbulence 

529  variances  respond  to  the  wave  motion  in  dramatic  ways.  The  horizontal  and  vertical  vari- 

530  ances  are  significantly  enhanced  by  the  motion  of  the  wave  surface  in  the  low-wind  cases. 

531  Even  though  the  turbulence  is  relatively  weak  the  turbulence  variances  are  large  near  the 

532  wave  surface  due  to  wave  pumping. 

533  6.  Summary 

534  A  large-eddy  simulation  (LES)  model  for  the  marine  atmospheric  planetary  boundary 

535  layer  (PBL)  is  coupled  to  a  3D  time-dependent  surface  gravity  wavefield.  A  coordinate 

536  transform  from  physical  to  computational  space  is  used  that  accounts  for  vertical  movement 

537  of  the  mesh  in  physical  space.  We  use  the  geometric  conservation  law  (GCL)  (Thomas  and 

538  Lombard  1979)  to  solve  for  the  grid  speeds  that  enter  into  the  advection  of  momentum 

539  and  scalars.  The  algorithm  is  used  to  carry  out  a  series  of  simulations  over  a  broadband 

540  moving  wavefield  that  conforms  to  a  Pierson-Moskowitz  wave  spectrum.  The  wave  age 

541  Cp/Uio  =  [1.5, 4.8]  varies  from  near  wind-wave  equilibrium  to  a  low-wind  swell  dominated 

542  regime.  In  the  low  wind  cases  we  find  features  similar  to  previous  observational  and  modeling 

543  investigations:  the  surface  layer  winds  show  clear  departures  from  rough  wall  law-of-the-wall 

544  scaling.  The  coherence  and  magnitude  of  the  pressure  field  p'/p  depends  critically  on  the 

545  motion  of  the  underlying  wavefield  and  the  turbulence  level. 
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APPENDIX 


557 

558  Appendix  Computation  of  surface  fluxes 

559  At  each  (x,y)  or  (£,77)  surface  point,  surface  fluxes  are  computed  in  a  local  wave-fitted 

560  Cartesian  coordinate  system  xd  assuming  law-of-the-wall  formulas.  These  fluxes  are  then 

561  converted  into  the  LES  Cartesian  system  x  so  that  they  can  be  used  as  SGS  surface  fluxes 

562  in  (9).  The  sequence  of  the  steps  is: 

563  1.  Build  the  base  vectors  parallel,  perpendicular,  and  normal  to  the  wave  surface  that 

564  define  the  local  xf  coordinate  system.  Compute  the  surface  wind 

565  2.  Build  the  matrix  of  direction  cosines  a  that  transforms  x!  x 

566  3.  Build  the  surface  momentum  and  temperature  fluxes  in  the  x!  system  based  on  Monin- 

567  Obukhov  (MO)  similarity  theory 

568  4.  Transform  the  momentum  fluxes  in  the  xd  coordinate  system  into  the  LES  coordinate 

569  system  x 


570  a.  Surface  definition  and  surface  wind 


571  Given  the  (x,  y ,  z)  surface  coordinates  the  vectors  aligned  with  the  surface  (£,  77)  gridlines 

572  for  (  =  0  are 

/\z  /\z 

573  t\  =  i  +  — — k  (constant  77)  ,  to  =  j  +  — — k  (constant  £)  ,  (Al) 

Ax  Ay 

574  where  (Ax,  Ay,  A z)  are  displacements  taken  along  the  surface  gridlines.  The  surface  normal 


fl  —  t\  X  t^  , 


577  and  the  unit  base  vectors  (i\  —  L/|fi  |,  h  —  h/lk],  n  —  n/|n|)  are  then  easily  computed. 
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578  In  the  simulation  code,  the  components  of  the  Cartesian  wind  vector  at  the  first  grid 

579  point  off  the  surface  are  u  =  ui  +  uj  +  wk,  and  in  the  situation  of  moving  water  waves  the 

580  velocity  components  of  the  surface  are  uQ  =  uQ i  +  vaj  +  «yk.  Following  standard  practice, 

581  we  assume  that  the  surface  stress  acts  in  a  direction  parallel  to  the  surface  and  depends  on 

582  the  relative  motion  between  the  wind  and  water  (e.g.,  Li  1995).  The  relative  wind  vector 

583  us  =  Su  =  u  —  ua  is  projected  onto  the  surface,  i.e., 

584  us  =  (Su  ■  ti  )t\  +  (Su  ■  t2)t2  =  u^ti  +  u^t2 ,  (A3) 

585  where  (ui1\ui2'))  are  the  relative  surface  wind  components  aligned  with  the  vectors  (t\ ,  t2), 

586  respectively. 


587  b.  Bulk  transfer  coefficients  Cd  and  Ch 


588  To  compute  the  magnitude  of  the  total  surface  stress  rG  we  invoke  MO  similarity  theory 

589  and  the  bulk  aerodynamic  formulas  for  momentum  and  heat: 

590  t0  =  —ul  =  —  Cd\us\2  (A4a) 

591  Q*  hAltisI  (0 s  @surf^)  •  (A4b) 


592  where  (Cd,  Ch)  are  bulk  transfer  coefficients,  w*  is  the  friction  velocity,  9S  is  the  potential 

593  temperature  at  the  first  grid  point  off  the  surface,  i.e.,  at  the  same  location  as  |us|,  0surf  is 

594  the  temperature  of  the  surface,  and  Q*  is  the  surface  heat  (temperature)  flux.  Under  the 

595  assumption  that  the  surface  shear  stress  points  in  the  direction  of  the  unit  base  vector  us: 


T  —  TaUs  —  CJ fj  I  Us 


(1)  (2) 
v>s  a  Ul  . 
rti  +  3 - rt2 


\us 


597  or 


598  T  =  -Cd\us\  (u^Hi  +  ufH2) 

599  For  a  flat  surface  with  =  0  then  (A6)  becomes 

600  T  =  —  Crf|«s|  (u^l  +  t42)j) 


(A5) 


(A6) 


(AT) 
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601  which  is  the  variation  in  flat  LES  codes.  In  the  xf  coordinate  system,  the  surface  momentum 

602  flux  r7  is  a  symmetric  matrix  with  two  non- zero  components 


603 


Tij  —  T° 


0  0  1 
0  0  0 
1  0  0 


(A8) 


604  Rules  need  to  be  formulated  for  the  transfer  coefficients.  In  the  absence  of  detailed  infor- 

605  mation  we  simply  use  MO  similarity  arguments  as  for  a  flat  surface.  Given  the  expressions 


606 

|«s| 

— 

■u* 

hi  ( Zm] 

|  +  V’m  ( y ) 

hi 

\ZoJ 

jL/ 

9,-9 


surf 


Q* 


(A9a) 

(A9b) 


60s  where  the  von  Karman  constant  is  n,  the  surface  roughness  is  zQ,  the  normal  distance  from 

609  the  surface  to  the  first  gridpoint  is  zs,  and  the  MO  length  is  L.  The  functions  ifh  are 

610  the  traditional  MO  similarity  functions  (e.g.,  Large  et  al.  1994).  We  use  iteration  to  solve 

611  (A9a,  A9b)  for  u*  and  9surf  when  Q*  is  given. 


612  c.  Matrix  of  direction  cosines  [a] 

613  The  matrix  of  direction  cosines  connects  the  xt  and  x  coordinate  systems: 

614  xf  —  ax  where  a  =  [us , f, n]  .  (A10) 

615  (us,f,n)  are  unit  base  vectors  with  f  =  h  x  us.  Transformation  of  stress  between  two 

616  Cartesian  coordinate  systems  obeys  the  transformation  rules  for  a  second-order  tensor  {e.g., 

617  see  p.  147,  Goldstein  1950) 

618  r  =  araT  or  aTra  =  r.  (All) 

619  Expanding  the  matrix  multiplications  in  (All)  and  utilizing  the  properties  of  at]  the  transfor- 

620  mation  of  stress  from  the  local  wave-fitted  Cartesian  coordinate  system  to  the  LES  coordinate 
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621  system  is 


622 


r  =  r0bij 


To 


2auCl3i  ®11®32  +  ®31®12  Olla33  +  a13°31 

«12«31  +  Q'll®32  2a\2a32  a  12«33  +  a13°32 

®11®33  +  ®13a31  «13a32  +  «12«33  2013033 


(A12) 


623  The  right  hand  side  of  (A12)  is  a  full  symmetric  matrix.  If  the  surface  is  2-D,  i.e.,  no 

624  variation  in  y  (or  rj),  the  direction  cosine  of  the  j-component  of  the  surface  normal  vector 

625  032  =  0.  Then  (A12)  reduces  to  a  similar  transformation  matrix  as  described  by  Henn  and 
Sykes  (1999)  (see  their  equation  (16)). 


626 


627 


628  REFERENCES 

629  Anderson,  D.  A.,  J.  C.  Tannehill,  and  R.  H.  Pletcher,  1984:  Computational  Fluid  Mechanics 

630  and  Heat  Transfer.  McGraw-Hill,  599  pp. 

631  Anderson,  W.  and  C.  Meneveau,  2010:  A  large-eddy  simualtion  model  for  boundary-layer 

632  flow  over  surfaces  with  horizontally  resolved  and  vertically  unresolved  roughness  elements. 

633  Boundary- Layer  Meteorol.,  137,  397-415. 

634  Armheld,  S.  and  R.  Street,  2002:  An  analysis  and  comparison  of  the  time  accuracy  of 

635  fractional-step  methods  for  the  Navier-Stokes  equations  on  staggered  grids.  International 

636  Journal  for  Numerical  Methods  In  Fluids ,  38,  255-282. 

637  Ayotte,  K.,  et  ah,  1995:  An  evaluation  of  neutral  and  convective  planetary  boundary  layer 

638  parameterizations  relative  to  large  eddy  simulation.  Boundary- Layer  Meteorol.,  79,  131— 

639  175. 

640  Ayotte,  K.  W.  and  D.  E.  Hughes,  2004:  Observations  of  boundary-layer  wind-tunnel  flow 

641  over  isolated  ridges  of  varying  steepness  and  roughness.  Boundary- Layer  Meteorol.,  112, 

642  525-556. 

643  Banner,  M.  L.,  1990:  The  influence  of  wave  breaking  on  the  surface  pressure  distribution  in 

644  wind-wave  interaction.  J.  Fluid  Mech.,  211,  463-495. 

645  Banner,  M.  L.  and  W.  K.  Melville,  1976:  On  the  separation  of  airflow  over  water  waves.  J. 

646  Fluid  Mech.,  77,  825-842. 

647  Beare,  R.  J.,  et  al.,  2006:  An  intercomparison  of  large-eddy  simulations  of  the  stable  bound- 

648  ary  layer.  Boundary-Layer  Meteorol,  118,  242-272. 


29 


649 


Belcher,  S.  E.  and  J.  C.  R.  Hunt,  1998:  Turbulent  flow  over  hills  and  waves.  Annual  Review 


650  of  Fluid  Mechanics ,  30,  507-538. 

651  Belcher,  S.  E.,  et  al.,  2012:  A  global  perspective  on  Langmuir  turbulence  in  the  ocean  surface 

652  boundary  layer.  Geophys.  Res.  Lett.,  39,  L18605. 

653  Black,  P.,  et  ah,  2007:  Air-sea  exchange  in  hurricanes:  Synthesis  of  observations  from  the 

654  Coupled  Boundary  Layers  Air-Sea  Transfer  experiment.  Bull.  Amer.  Meteorol.  Soc.,  88, 

655  357-374. 

656  Brown,  A.  R.,  J.  M.  Hobson,  and  N.  Wood,  2001:  Large  eddy  simulation  of  neutral  turbulent 

657  flow  over  rough  sinusoidal  ridges.  Boundary- Layer  Meteorol.,  98,  411-441. 

658  Chalikov,  D.,  1998:  Interactive  modeling  of  surface  waves  and  boundary  layer.  Ocean  Wave 

659  Measurement  and  Analysis:  Proceedings  of  the  Third  International  Symposium  Waves  97, 

660  B.  L.  Edge,  J.  M.  Hemsley,  and  Y.  Goda,  Eds.,  American  Society  of  Civil  Engineers, 

661  1525-1539. 

662  Chalikov,  D.  and  S.  Rainchik,  2011:  Coupled  numerical  modelling  of  winds  and  waves  and 

663  the  theory  of  the  wave  boundary  layer.  Boundary- Layer  Meteorol.,  138,  1-41. 

664  Chalikov,  D.  V.,  1978:  The  numerical  simulation  of  wind-wave  interaction.  J.  Fluid  Mech ., 

665  87,  561-582. 

eee  Chen,  S.  S.,  J.  F.  Price,  W.  Zhao,  M.  A.  Donelan,  and  E.  J.  Walsh,  2007:  The  CBLAST- 
667  Hurricane  program  and  the  next-generation  fully  coupled  atmosphere-wave-ocean  models 
ees  for  hurricane  research  and  prediction.  Bull.  Amer.  Meteorol.  Soc.,  88,  311-317. 

669  Deardorff,  J.  W.,  1972:  Three-dimensional  numerical  modeling  of  the  planetary  boundary 

670  layer.  Workshop  on  Micrometeorology,  D.  A.  Haugen,  Eel.,  American  Meteorological  Soci- 

671  ety,  271-311. 


30 


672 


Demirdzic,  I.  and  M.  Peric,  1990:  Finite  volume  method  for  prediction  of  fluid  flow  in 

673  arbitrarily  shaped  domains  with  moving  boundaries.  International  Journal  for  Numerical 

674  Methods  in  Fluids,  10,  771-790. 

675  Donelan,  M.  A.,  1998:  Air-water  exchange  processes.  Physical  Processes  in  Lakes  and  Oceans, 

676  J.  Imberger,  Ed.,  American  Geophysical  Union,  Coastal  and  Estuarine  Studies,  Vol.  54, 

677  19—36. 

678  Donelan,  M.  A.,  J.  Hamilton,  and  W.  H.  Hui,  1985:  Directional  spectra  of  wind-generated 

679  waves.  Phil.  Trans.  Roy.  Soc.  London,  Ser.  A,  315,  509-562. 

680  Ferziger,  J.  H.  and  M.  Peric,  2002:  Computational  Methods  for  Fluid  Dynamics.  3d  ed., 

681  Springer- Verlag. 

682  Gent,  P.  R.,  1977:  A  numerical  model  of  the  air  flow  above  water  waves.  Part  2.  J.  Fluid 

683  Mech.,  82,  349-369. 

684  Gent,  P.  R.  and  P.  A.  Taylor,  1976:  A  numerical  model  of  the  air  flow  above  water  waves, 

ess  J.  Fluid  Mech.,  77,  105-128. 

ese  Ghosal,  S.  and  P.  Moin,  1995:  The  basic  equations  for  the  large  eddy  simulation  of  turbulent 
687  flows  in  complex  geometry.  J.  Comp.  Phys.,  118,  24-37. 

ess  Goldstein,  H.,  1950:  Classical  Mechanics.  Addison- Wesley  Publishing  Company,  399  pp. 

ess  Gong,  W.,  P.  A.  Taylor,  and  A.  Dornbrack,  1996:  Turbulent  boundary-layer  flow  over  fixed 
ego  aerodynamically  rough  two-dimensional  sinusoidal  waves.  J.  Fluid  Mech.,  312,  1-37. 

691  Grachev,  A.  A.  and  C.  W.  Fairall,  2001:  Upward  momentum  transfer  in  the  marine  boundary 

692  layer.  J.  Phys.  Oceanogr.,  31,  1698-1711. 

693  Hanley,  K.  E.  and  S.  E.  Belcher,  2008:  Wave-Driven  wind  jets  in  the  marine  atmospheric 

694  boundary  layer.  J.  Atmos.  Sci.,  65,  2646-2660. 


31 


695  Hanley,  K.  E.,  S.  E.  Belcher,  and  P.  P.  Sullivan,  2010:  A  global  climatology  of  wind-wave 

696  interaction.  J.  Phys.  Oceanogr .,  40,  1263-1282. 

697  Harris,  J.  A.,  S.  E.  Belcher,  and  R.  L.  Street,  1996:  Linear  dynamics  of  wind  waves  in 

698  coupled  turbulent  air-water  flow.  Part  2.  Numerical  model.  J.  Fluid  Mech.,  308,  219-254. 

699  Henn,  D.  S.  and  R.  I.  Sykes,  1999:  Large-eddy  simulations  of  flow  over  wavy  surfaces.  J. 

700  Fluid  Mech.,  383,  75-112. 

701  Kornen,  G.  J.,  L.  Cavaleri,  M.  Donelan,  K.  Hasselmann,  S.  Hasselmann,  and  P.  A.  E.  M. 

702  Janssen,  1994:  Dynamics  and  Modelling  of  Ocean  Waves.  Cambridge  LIniversity  Press, 

703  532  pp. 

704  Large,  W.  G.,  J.  C.  McWilliams,  and  S.  C.  Doney,  1994:  Oceanic  vertical  mixing:  A  review 

705  and  a  model  with  a  nonlocal  boundary  layer  parameterization.  Rev.  Geophys.,  32,  363-403. 

706  Le,  H.  and  P.  Moin,  1991:  An  improvement  of  fractional  step  methods  for  the  incompressible 

707  Navier-Stokes  equations.  J.  Comp.  Phys.,  92,  369-379. 

708  Li,  P.  Y.,  1995:  A  numerical  study  on  energy  transfer  between  turbulent  air  flow  and  finite 

709  amplitude  water  waves.  Ph.D.  thesis,  York  LIniversity,  182  pp.,  North  York,  Ontario. 

710  Lighthill,  J.,  1978:  Waves  in  Fluids.  Cambridge  University  Press,  504  pp. 

711  Lundquist,  K.  A.,  F.  K.  Chow,  and  J.  K.  Lundquist,  2010:  An  immersed  boundary  method 

712  for  the  Weather  Research  and  Forecasting  model.  Mon.  Wea.  Rev.,  138,  796-817. 

713  Makin,  V.  K.,  2008:  On  the  possible  impact  of  a  following-swell  on  the  atmospheric  boundary 

714  layer.  Boundary- Layer  Meteorol.,  129,  469-478. 

715  Makin,  V.  K.,  V.  N.  Kudryavtsev,  and  C.  Mastenbroek,  1995:  Drag  of  the  sea  surface. 

716  Boundary- Layer  Meteorol.,  73,  159-182. 


32 


717 


McWilliams,  J.  C.,  E.  Huckle,  J.-H.  Liang,  and  P.  P.  Sullivan,  2012:  The  wavy  Ekman  layer: 

718  Langmuir  circulations,  breakers  and  Reynolds  stress.  J.  Phys.  Oceanogr.,  42,  1793-1816. 

719  McWilliams,  J.  C.,  C.-H.  Moeng,  and  P.  P.  Sullivan,  1999:  Turbulent  fluxes  and  coherent 

720  structures  in  marine  boundary  layers:  Investigations  by  large-eddy  simulation.  Air-Sea 

721  Exchange:  Physics,  Chemistry,  Dynamics,  and  Statistics ,  G.  Geernaert,  Ed.,  Kluwer, 

722  507-538. 

723  Meirink,  J.  F.  and  V.  K.  Makin,  2001:  Modelling  low-Reynolds-number  effects  in  the  turbu- 

724  lent  air  flow  over  water  waves.  J.  Fluid  Mech.,  415,  155-174. 

725  Mittal,  R.  and  G.  Iaccarino,  2005:  Immersed  boundary  methods.  Annual  Review  of  Fluid 

726  Mechanics,  37,  239-261. 

727  Moeng,  C.-H.,  1984:  A  large-eddy  simulation  model  for  the  study  of  planetary  boundary- 

728  layer  turbulence.  J.  Atmos.  Sci.,  41,  2052-2062. 

729  Moeng,  C.-H.,  P.  P.  Sullivan,  M.  F.  Khairoutdinov,  and  D.  A.  Randall,  2010:  A  mixed 

730  scheme  for  subgrid-scale  fluxes  in  cloud  resolving  models.  J.  Atmos.  Sci.,  67,  3692-3705. 

731  Moeng,  C.  H.  and  J.  C.  Wyngaard,  1988:  Spectral  analysis  of  large-eddy  simulations  of  the 

732  convective  boundary  layer.  J.  Atmos.  Sci.,  45,  3573-3587. 

733  Nakayama,  A.,  K.  Hori,  and  R.  L.  Street,  2004:  Filtering  and  LES  of  flow  over  irregular 

734  rough  boundary.  Center  for  Turbulence  Research,  Proceedings  of  the  Summer  Program, 

735  Stanford  University/NASA  Ames  Research  Center,  145-156. 

736  Nakayama,  A.,  K.  Sakio,  Y.  Kitano,  and  S.  Yokojima,  2009:  Direct  numerical  simulation 

737  and  filtering  of  turbulent  flows  over  model  rough  surfaces.  Memoirs  of  the  Graduate  School 

738  of  Engineering  Kobe  University,  1,  9-41. 

739  Nilsson,  E.  O.,  A.  Rutgersson,  A.  Smedman,  and  P.  Sullivan,  2012:  Convective  boundary 


33 


740  layer  structure  in  the  presence  of  wind-following  swell.  Quart.  J.  Roy.  Meteorol.  Soc .,  138, 

741  1476-1489. 

742  Phillips,  O.  M.,  1977:  Dynamics  of  the  Upper  Ocean.  Cambridge  University  Press,  336  pp. 

743  Rhie,  C.  M.  and  W.  L.  Chow,  1983:  A  numerical  study  of  the  turbulent  flow  past  an  isolated 

744  airfoil  with  trailing  edge  separation.  AIAA  J.,  21,  1525-1532. 

745  Romero,  L.  and  W.  K.  Melville,  2010:  Airborne  observations  of  fetch-limited  waves  in  the 

746  Gulf  of  Tehuantepec.  J.  Phys.  Oceanogr.,  40,  441-465. 

747  Sajjadi,  S.  G.,  J.  C.  R.  Hunt,  and  F.  Drullion,  2013:  Asymptotic  multi-layer  analysis  of  wind 

748  over  unsteady  monochromatic  surface  waves.  J.  of  Engineering  Mathematics ,  in  press. 

749  Smedman,  A.,  U.  Hogstrom,  H.  Bergstrom,  and  A.  Rutgersson,  1999:  A  case  study  of  air-sea 

750  interaction  during  swell  conditions.  J.  Geophys.  Res.,  104,  25  833-25  851. 

751  Spalart,  P.  R.,  R.  D.  Moser,  and  M.  M.  Rogers,  1991:  Spectral  methods  for  the  Navier-Stokes 

752  equations  with  one  infinite  and  two  periodic  directions.  J.  Comp.  Phys.,  96,  297. 

753  Sullivan,  P.  P.,  J.  B.  Edson,  T.  Hristov,  and  J.  C.  McWilliams,  2008:  Large  eddy  simulations 

754  and  observations  of  atmospheric  marine  boundary  layers  above  non-equilibrium  surface 

755  waves.  J.  Atmos.  SO.,  65,  1225-1245. 

756  Sullivan,  P.  P.,  T.  W.  Horst,  D.  H.  Lenschow,  C.-H.  Moeng,  and  J.  C.  Weil,  2003:  Structure 

757  of  subfilter-scale  fluxes  in  the  atmospheric  surface  layer  with  application  to  large-eddy 

758  simulation  modeling.  J.  Fluid  Mech.,  482,  101-139. 

759  Sullivan,  P.  P.  and  J.  C.  McWilliams,  2002:  Turbulent  flow  over  water  waves  in  the  presence 

760  of  stratification.  Phys.  Fluids,  14,  1182-1195. 

761  Sullivan,  P.  P.  and  J.  C.  McWilliams,  2010:  Dynamics  of  winds  and  currents  coupled  to 

762  surface  waves.  Annual  Review  of  Fluid  Mechanics,  42,  19-42. 


34 


763  Sullivan,  P.  P.,  J.  C.  McWilliams,  and  C.-H.  Moeng,  1996:  A  grid  nesting  method  for 

764  large-eddy  simulation  of  planetary  boundary  layer  flows.  Boundary- Layer  Meteorol.,  80, 

765  167-202. 

766  Sullivan,  P.  P.,  J.  C.  McWilliams,  and  C.-H.  Moeng,  2000:  Simulation  of  turbulent  flow  over 

767  idealized  water  waves.  J.  Fluid  Mech.,  404,  47-85. 

768  Sullivan,  P.  P.  and  E.  G.  Patton,  2011:  The  effect  of  mesh  resolution  on  convective  boundary- 

769  layer  statistics  and  structures  generated  by  large-eddy  simulation.  J.  Atmos.  Sci.,  68, 

770  2395-2415. 

771  Sullivan,  P.  P.,  E.  G.  Patton,  and  K.  W.  Ayotte,  2010:  Turbulent  flow  over  and  around 

772  sinusoidal  bumps,  hills,  gaps  and  craters  derived  from  large  eddy  simulations.  19th  Amer. 

773  Meteorol.  Soc.  Symp.  on  Boundary  Layer  and  Turbulence,  Keystone,  CO. 

774  Sullivan,  P.  P.,  L.  Romero,  J.  C.  McWilliams,  and  W.  K.  Melville,  2012:  Transient  evolution 

775  of  Langmuir  turbulence  in  ocean  boundary  layers  driven  by  hurricane  winds  and  waves. 

776  J.  Phys.  Oceanogr .,  42,  1959-1980. 

777  Taylor,  P.  A.,  1998:  Turbulent  boundary-layer  flow  over  low  and  moderate  slope  hills.  Journal 

778  of  Wind  Engr.  and  Industrial  Aerodyanmics ,  74-76,  25-47. 

779  Thomas,  P.  D.  and  C.  K.  Lombard,  1979:  Geometric  conservation  law  and  its  application  to 

780  flow  computations  on  moving  grids.  AIAA  Journal,  17,  1030-1037. 

781  Wilson,  J.  D.,  2002:  Representing  drag  on  unresolved  terrain  as  a  distributed  momentum 

782  sink.  J.  Atmos.  Sci.,  59,  1629-1637. 

783  Wyngaard,  J.  C.,  1998:  Boundary-Layer  modeling:  History,  philosophy,  and  sociology.  Clear 

784  and  Cloudy  Boundary  Layers,  A.  A.  M.  Holtslag  and  P.  G.  Duynkerke,  Eds.,  Royal  Nether- 

785  lands  Academy  of  Arts  and  Sciences. 


35 


786 


Wyngaard,  J.  C.,  2004:  Toward  numerical  modeling  in  the  Terra  Incognita.  J.  Atmos.  Sci ., 


ysr  61,  1816-1826. 

788  Zang,  Y.,  R.  L.  Street,  and  J.  R.  Koseff,  1994:  A  non-staggered  grid,  fractional  step  method 

789  for  time-dependent  incompressible  Navier-Stokes  equations  in  curvilinear  coordinates.  J. 

790  Comp.  Phys.,  114,  18-33. 

791  Zheng,  Z.  and  L.  Petzold,  2006:  Runge-Kutta-Chebyshev  projection  method.  J.  Comp. 

792  Phys.,  219,  976-991. 


36 


793 


List  of  Tables 


794  1  Simulation  properties 


38 


37 


Table  1:  Simulation  properties 


Run 

Ug  (m  s  x) 

cp/u10 

■u*  (m  s  *) 

A 

5 

4.8 

0.124 

B 

7.5 

3.4 

0.187 

C 

10 

2.8 

0.228 

D 

15 

1.9 

0.338 

E 

20 

1.5 

0.452 
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3  A  snapshot  of  the  wavefield  height  h(x,  y,  t )  that  is  imposed  at  the  bottom 
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5  An  instantaneous  x  —  z  slice  of  the  3D  time  varying  computational  mesh  in 
the  lowest  portion  of  the  PBL.  The  (£,  rj)  gridlines  become  level  surfaces  at 
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Snapshot  of  static  pressure  fluctuations  p' / p  in  an  x  —  z  plane  near  the  water 
surface.  The  upper  panel  is  a  swell  dominated  regime  with  wave  age  ~  4.8 
while  the  lower  panel  is  a  case  near  wind-wave  equilibrium  with  wave  age 
~  1.4.  The  wave  spectrum  is  a  Pierson-Moskowitz  spectrum.  Notice  the 
coherence  between  the  wave  held  and  the  pressure  fluctuations  in  the  case 
with  swell.  The  color  bar  is  in  units  of  m  s~2  and  the  range  is  different 
between  the  two  cases. 

Snapshot  of  resolved  vertical  velocity  fluctuations  liZ/w*  in  a  wave  following 
x  —  y  plane  near  the  water  surface  (  =  2.5  m.  The  left  panel  is  a  swell 
dominated  regime  with  wave  age  ~  4.8  while  the  right  panel  is  a  case  near 
wind-wave  equilibrium  with  wave  age  ~  1.4.  The  wave  spectrum  at  the  bot¬ 
tom  of  the  PBLs  is  the  same.  Notice  the  range  of  the  color  bar  is  different 
between  the  two  cases.  The  (normalized)  fluctuations  in  the  wind-wave  equi¬ 
librium  case  are  smaller. 

Vertical  profiles  of  wind  speed  (left  panel)  and  turbulence  variances  (right 
panels)  for  different  values  of  wave  age  Cp/U\q.  Friction  velocity  w*  is  used 
for  normalization.  The  dashed  black  line  is  the  rough  wall  formula  U/u*  — 
In (z/z0)/k,  where  k  =  0.4.  Temporal  and  spatial  averaging  is  used  to  make 
the  statistics.  The  spatial  averaging  is  over  £  —  77  planes,  i.e.,  at  constant  £. 


Figure  1:  A  two-dimensional  sketch  illustrating  the  layout  of  the  cell  centered  Cartesian 
velocity  components  ( u,w ),  along  with  the  pressure,  potential  temperature,  and  SGS  energy 
(p*,  9 ,  e).  The  contravariant  flux  velocity  W  is  located  at  a  cell  face  and  is  oriented  perpen¬ 
dicular  to  a  (  =  constant  surface.  The  contravariant  flux  velocity  U  is  oriented  perpendicular 
to  a  £  =  constant  surface  and  is  also  located  at  a  cell  center,  see  Section  3a.  The  horizontal 
velocity  components  (v,  V),  which  are  not  shown,  point  into  the  page  at  the  cell  center. 
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Figure  2:  A  two-dimensional  sketch  illustrating  the  layout  of  the  cell  centered  pressure  p* 
near  the  wavy  boundary  (  =  0.  The  pairs  of  indices  [(i  —  1,  0), (i  +  1,  2)]  are  the  discrete 
locations  of  the  cell  centers  in  the  (£,  ()  directions,  respectively.  Ghost  point  pressures  are 
located  along  lines  of  constant  (  =  —A(/2,  i.e.,  along  lines  (i,0),i  =  1  ,NX. 
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Figure  3:  A  snapshot  of  the  waveheld  height  h(x,y,t)  that  is  imposed  at  the  bottom  of  the 
LES  code.  The  horizontal  grid  spacing  matches  the  LES,  i.e.,  Ax  =  Ay  =  2.93  m.  The 
color  bar  is  in  units  of  meters. 
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Figure  4:  Normalized  1-D  power  spectra  of  the  wave  height  as  functions  of  the  horizontal 
wavenumbers  E(kx )  and  E(ky )  normalized  by  the  total  variance  ( h 2). 
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Figure  5:  An  instantaneous  x  —  z  slice  of  the  3D  time  varying  computational  mesh  in  the 
lowest  portion  of  the  PBL.  The  (£,77)  gridlines  become  level  surfaces  at  about  100  m  above 
the  water.  Only  a  fraction  of  the  grid  is  displayed. 
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Figure  6:  Snapshot  of  static  pressure  fluctuations  p' / p  in  an  x  —  z  plane  near  the  water 
surface.  The  upper  panel  is  a  swell  dominated  regime  with  wave  age  ~  4.8  while  the  lower 
panel  is  a  case  near  wind-wave  equilibrium  with  wave  age  ~  1.4.  The  wave  spectrum  is  a 
Pierson-Moskowitz  spectrum.  Notice  the  coherence  between  the  wave  field  and  the  pressure 
fluctuations  in  the  case  with  swell.  The  color  bar  is  in  units  of  m  s~2  and  the  range  is 
different  between  the  two  cases. 
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Figure  7:  Snapshot  of  resolved  vertical  velocity  fluctuations  w'/u*  in  a  wave  following  x  —  y 
plane  near  the  water  surface  (  =  2.5  m.  The  left  panel  is  a  swell  dominated  regime  with 
wave  age  ~  4.8  while  the  right  panel  is  a  case  near  wind-wave  equilibrium  with  wave  age 
~  1.4.  The  wave  spectrum  at  the  bottom  of  the  PBLs  is  the  same.  Notice  the  range  of  the 
color  bar  is  different  between  the  two  cases.  The  (normalized)  fluctuations  in  the  wind-wave 
equilibrium  case  are  smaller. 
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Figure  8:  Vertical  profiles  of  wind  speed  (left  panel)  and  turbulence  variances  (right  panels) 
for  different  values  of  wave  age  Cp/Uiq.  Friction  velocity  w*  is  used  for  normalization.  The 
dashed  black  line  is  the  rough  wall  formula  U/u*  =  In (z/z0)/k,,  where  k  =  0.4.  Temporal 
and  spatial  averaging  is  used  to  make  the  statistics.  The  spatial  averaging  is  over  £  —  r) 
planes,  i.e.,  at  constant  (. 
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4B.4  A  LARGE  EDDY  SIMULATION  MODEL  OF  HIGH  WIND  MARINE 

BOUNDARY  LAYERS  ABOVE  A  SPECTRUM  OF  RESOLVED  MOVING  WAVES 
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1.  INTRODUCTION 

An  improved  understanding  of  the  interaction  of 
winds,  waves,  and  currents  in  the  upper  ocean  at  scales 
of  individual  waves  and  wave  groups  is  needed  to  fur¬ 
ther  develop  the  next  generation  of  climate,  weather,  and 
wave  forecast  models.  For  example,  coupled  wind-wave- 
ocean  models  (Chen  et  al.,  2007;  Black  et  ah,  2007)  are 
viewed  as  critical  tools  for  accurate  prediction  of  tropi¬ 
cal  cyclone  intensity  and  track  forecasts,  but  these  mod¬ 
eling  systems  employ  a  suite  of  parameterizations  that 
are  largely  statistical  descriptions  of  the  wind-wave  in¬ 
teractions  that  generate  the  critical  momentum  and  scalar 
fluxes.  These  forecast  models  do  not  account  for  the  im¬ 
portant  phase  relationships  between  winds,  waves  and 
currents,  e.g.,  the  spatial  and  temporal  intermittency  of 
wave  breaking  that  occurs  in  moderate  to  high  winds 
(see  figure  1).  Also,  there  is  a  growing  appreciation 
that  wave-current  interactions  are  important  for  the  upper 
ocean  boundary  layer  (Sullivan  and  McWilliams,  2010), 
and  thus  for  climate  predictions,  and  that  remotely  gener¬ 
ated  swell  and  non-equilibrium  wave  states  can  play  an 
important  and  critical  role  in  the  surface-layer  dynam¬ 
ics  of  the  atmospheric  planetary  boundary  layer  (PBL) 
(Hanley  et  ah,  2010;  Sullivan  and  McWilliams,  2010). 

The  present  work  describes  our  recent  developments 
in  coupling  a  turbulence-resolving  large-eddy  simulation 
(LES)  of  the  atmospheric  PBL  to  a  three-dimensional 
time-dependent  resolved  surface  wavefield.  This  builds 
on  our  past  efforts  which  coupled  a  turbulent  boundary- 
layer  flow  with  a  single  monochromatic  wave  (Sullivan 
et  ah,  2000;  Sullivan  and  McWilliams,  2002;  Sullivan 
et  ah,  2008).  The  computational  method  described  here 
allows  for  nearly  arbitrary  3D  wavefields,  i.e.,  the  sea 
surface  elevation  h  =  h(x,y,t),  as  a  surface  boundary 
condition.  The  spatial  scales  of  the  resolved  turbulence 
and  waves  are  O(m)  up  to  the  scale  of  the  PBL  height. 
At  present,  the  waves  are  externally  imposed  based  on 
empirical  wave  spectra.  Ultimately  the  wavefields  will 

*  corresponding  author  address:  Peter  P.  Sullivan,  National  Center 
for  Atmospheric  Research,  P  O.  Box  3000,  Boulder,  CO  80307-3000; 
email :  pps@ucar.edu 


be  direct  observations  of  the  sea  surface  from  field  cam¬ 
paigns.  High  resolution  simulations  of  PBL  turbulence 
in  the  presence  of  surface  waves  has  the  potential  to  pro¬ 
vide  new  insight  into  the  dynamics  of  air-sea  coupling 
at  small  scales.  The  jump  in  massively  parallel  compu¬ 
tational  resources  facilitates  the  coupling  of  winds  and 
waves  in  turbulence  resolving  simulations. 

2.  LES  ALGORITHM  WITH  MOVING  WAVES 

An  LES  model  for  an  atmospheric  PBL  with  a  sta¬ 
tionary  undulating  lower  boundary  is  described  in  Sul¬ 
livan  et  ah  (2010).  They  outline  the  algorithm,  numer¬ 
ical  details,  code  parallelization,  and  present  results  for 
boundary-layer  flows  over  and  around  two-  and  three- 
dimensional  orography.  Here  we  build  on  those  develop¬ 
ments  but  focus  on  the  additional  complications  caused 
by  the  temporal  movement  of  the  3D  lower  surface.  In 
the  description  of  the  model  equations,  given  in  Sec¬ 
tion  2.2,  the  following  notation  is  used:  pu  denotes  the 
Cartesian  components  of  momentum,  0  is  virtual  poten¬ 
tial  temperature,  e  is  the  subgrid-scale  energy,  and  II  is 
the  pressure  variable.  Quantities  with  an  overbar  (  )  are 
interpreted  as  LES  spatially  filtered  variables. 

2.1  Coordinate  transformation 

We  follow  Sullivan  et  ah  (2010)  and  adapt  our 
LES  model  with  a  flat  bottom  to  the  situation  with  a 
three-dimensional  time-dependent  boundary  shape  h  = 
h(x,y,t)  by  applying  a  transformation  to  the  physical 
space  coordinates  (x,y,z)  that  maps  them  onto  compu¬ 
tational  coordinates  (^,q,Q.  The  computational  mesh 
in  physical  space  is  surface  following,  non-orthogonal, 
and  time  varying.  Also,  vertical  gridlines  are  held  fixed 
at  a  particular  (x,y)  location  on  the  surface  but  the  lines 
are  permitted  to  undergo  vertical  translation  as  a  func¬ 
tion  of  time  t,  i.e.,  vertical  gridlines  are  wave  follow¬ 
ing.  The  transformation  which  obeys  these  constraints 
and  maps  the  physical  domain  to  a  flat  computational  do¬ 
main  x  =>■  E,  is  the  rule: 


Figure  1:  Photograph  of  the  sea  surface  generated  by  winds  of  approximately  15  m  s-1  during  the  High  Resolution 
Air-Sea  Interaction  (Hi-Res)  field  campaign  carried  out  in  June  2010.  Notice  the  extensive  white  capping  generated 
by  large-scale  plunging  and  spilling  breakers  and,  in  the  foreground,  a  web  of  small-scale  breakers.  The  photograph 
is  taken  from  the  R/V  Floating  Instmment  Platform  (FLIP)  courtesy  of  Tihomir  Hristov.  For  a  description  of  Hi-Res 
see  http://airsea.ucsd.edu/hires/. 
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The  differential  metrics  dxt/d^j  and  d^j/dxj, 

which  are 

needed  in  formulating  the  LES  model,  are 

connected 
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rules  in  (1)  we  have  the  reduced  set  of  non-zero  metric 
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t* 
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where  J  is  the  Jacobian.  The  time  dependence  of  the 
mapping  appears  in  (2a)  where  Zt  =  dz/dt  is  the  grid 
speed,  i.e.,  the  vertical  velocity  of  individual  gridpoints. 


2.2  LES  equations  in  curvilinear  coordinates  with 
time  dependence 

The  LES  equations  in  curvilinear  coordinates  are  de¬ 
rived  in  a  straightforward  manner  by  applying  the  chain 
rule  for  differentiation  but  we  pay  close  attention  to  the 
transformation  of  the  material  derivative  D()/Dt  and  the 
mass  conservation  equation.  The  set  of  LES  equations  in 
computational  coordinates  under  the  transformation  (1) 
and  (2)  are: 
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The  equations  are  expressed  in  strong  conservation  form 
using  the  “contravariant  flux”  velocity 


(3a)  is  the  mass  conservation  (continuity)  equation,  (3c) 
is  the  momentum  transport  equation,  (3d)  is  the  scalar 
transport  equation,  (3e)  is  the  subgrid-scale  energy  trans¬ 
port  equation,  and  (3f)  is  the  pressure  Poisson  equa¬ 
tion.  The  right  hand  sides  of  (3)  model  physical  pro¬ 
cesses  in  the  marine  atmospheric  PBL,  e.g.,  pressure 
gradients,  Coriolis  rotation,  divergence  of  subgrid-scale 
fluxes,  buoyancy,  and  in  the  case  of  the  SGS  e  equation 
also  diffusion  and  dissipation. 

The  time  dependence  of  the  grid  modifies  the  LES 
equations:  the  Jacobian  appears  inside  the  time  tendency 
of  each  transport  equation  and  as  expected  advection 
contains  a  contribution  from  the  grid  movement,  i.e.,  the 
total  vertical  flux  of  variable  \|/  depends  on  the  difference 
between  the  physical  velocity  and  grid  speed  (W  —  zf)\|/. 
Also,  as  a  consequence  of  writing  the  equations  in  strong 
conservation  form  we  need  to  satisfy  (3b)  which  is  a  sim¬ 
plified  form  of  the  so-called  geometric  or  space  conser¬ 
vation  law  (GCL)  first  discussed  by  Thomas  and  Lom¬ 
bard  (1979).  In  our  derivation,  (3b)  follows  directly  from 
considering  the  unsteady  form  of  the  mass  conservation 
equation  written  in  differential  form.  Thomas  and  Lom¬ 
bard  (1979)  and  also  Demirdzic  and  Peric  (1990)  derive 
the  GCL  by  considering  an  integral  form  of  the  equa¬ 
tions  of  motion.  Inspection  of  (3)  shows  that  if  the  ve¬ 
locity  and  scalar  fields  are  set  to  constant  values  then 
the  left  hand  sides  of  (3c-e)  reduce  to  (3b).  Hence  our 
numerical  method  needs  to  satisfy  our  reduced  form  of 
the  GCL  discretely  in  order  to  prevent  artificial  sources 
and  sinks  from  developing  in  the  computational  domain. 
Thus,  the  mixed  pseudospectral  finite-differencing  spa¬ 
tial  differencing  evaluates  the  advective  terms  in  (3c-e) 
in  flux  form  and  not  the  rotational  form  used  by  Moeng 
(1984)  or  skew-symmetric  form  used  by  Sullivan  et  al. 
(2000). 

We  use  (3b)  with  our  Runge-Kutta  time  stepping 
scheme  in  the  following  way:  Given  the  wave  height  dis¬ 
tribution  at  any  future  timestep,  say  tn+l,  we  first  pick  the 
vertical  distribution  of  gridpoints  z(lj,T|,  based  on 

the  computational  domain  height  and  h(£,,r\,tn+l).  We 
then  insert  the  Runge-Kutta  time  stepping  discretization 
in  (3b)  and  solve  for  the  required  matching  grid  speed 
Zt .  Since  our  time  stepping  is  explicit  the  value  of  Zt  di¬ 
agnosed  applies  at  the  time  level  tn.  In  other  words,  if 
the  wavefield  is  imposed  knowledge  of  its  future  shape 
is  needed  to  diagnose  the  grid  speed  at  an  earlier  time. 
This  is  analogous  to  our  grid  nesting  scheme  where  the 
vertical  velocity  at  a  future  time  is  imposed  on  a  nested 
boundary  (Sullivan  et  al.,  1996). 


We  emphasize  that  the  upper  boundary  of  our  compu¬ 
tational  domain  is  far  from  the  lower  surface  and  asymp¬ 
totes  to  a  flat  level  surface.  As  a  result  Zt  varies  with 
distance  from  the  wave  surface.  This  is  different  than 
the  scheme  proposed  by  Chalikov  (1998)  where  the  grid 
speed  is  equal  to  the  vertical  motion  of  the  wavefield  for 
all  z.  This  implies  that  the  computational  mesh  at  any  C 
including  the  upper  boundary,  mimics  the  shape  of  the 
lower  surface.  At  the  lower  boundary  z  =  h,  the  time  rate 
of  change  of  the  wavy  surface  is 

Zt  =  h,  =  w0  -  htu0  -  h^Vo  (5) 

where  (m0,v0,w0)  are  the  orbital  motions  of  the  wave- 
field.  The  kinematical  boundary  condition  is  then  W  — 
Zt  =  0  so  that  there  is  no  flow  through  the  boundary.  No¬ 
tice  the  definition  of  W(z  =  h)  given  by  (4)  is  consistent 
with  (5).  The  surface  grid  movement  also  leads  to  sub¬ 
tleties  in  formulating  the  pressure  Poisson  equation  and 
its  boundary  conditions. 

2.3  Wavefield  prescription 

To  complete  our  marine  PBL  LES  we  need  to  pre¬ 
scribe  the  surface  wavefield.  For  the  present  computa¬ 
tions  we  use  empirical  two-dimensional  statistical  wave 
spectra  (e.g.,  see  Komen  et  al.,  1994,  p.  187) 

E(k,$)  =  S{k)  D(k,$)  (6) 

where  the  amplitude  S(k)  and  directional  Dik.  ())  spectra 
depend  on  the  wavenumber  k  —  |k|  =  kx\  +  ky\ \ ,  wave 
direction  (])  and  surface  wind  speed  U\q.  For  the  ampli¬ 
tude  spectrum  we  choose  the  classic  Pierson-Moskowitz 
shape  (Pierson  and  Moskowitz,  1964;  Alves  et  al.,  2003) 
while  a  simple  directional  spectrum  is  designed  to  em¬ 
phasize  long-crested  waves  with  their  spanwise  axis  ori¬ 
ented  perpendicular  to  the  wind  direction  D(k,<\>)  =  (k  • 
U)".  In  physical  space,  the  synthetic  wavefield  h(x,y,t) 
is  constructed  from  a  sum  of  plane  waves 

h(x,t)  =^/z(k)exp[t(k-x  — COf)]  (7) 

k 

with  the  wave  amplitudes  /z(k)  picked  to  match  E(k,<\>). 
The  phases  are  chosen  from  a  Gaussian  distribution.  In¬ 
side  the  LES,  the  wavefield  is  advanced  in  time  using 
the  linear  dispersion  relation  or  =  gk  once  the  initial 
distribution  of  amplitudes  is  specified.  (7)  is  efficiently 
evaluated  using  2D  Fast  Fourier  Transforms.  Figure  2 
shows  a  typical  instantaneous  3D  wavefield  that  is  in¬ 
put  at  the  bottom  of  the  LES.  Future  computations  will 
use  the  detailed  phase-resolved  wave  measurements  col¬ 
lected  in  the  Hi-Res  field  campaign. 


3.  SIMULATIONS 


A  series  of  simulations  with  varying  geostrophic 
wind  Ug  =  (5,7.5,10,15,20)  m  s”1  are  carried  out 
for  a  neutrally-stratified  marine  PBL  in  a  domain 
(Xl.Yl,Zl)  =  (1200,1200,800)  m  using  (Nx,Ny,Nz)  = 
(512,512, 128)  gridpoints;  thus  the  horizontal  grid  spac¬ 
ing  Ax  =  Ay  =  2.34  m  and  the  first  vertical  level  is  1 
m  above  the  water.  A  slice  of  the  computational  mesh 
is  given  in  figure  3.  The  initial  temperature  sound¬ 
ing  0  =  300  K  up  to  the  inversion  height  Zi  =  400  m, 
beyond  this  height  0  increases  linearly  at  3  x  10^3  K 
m_1.  The  surface  heating  Q*  =  0,  the  surface  rough¬ 
ness  z0  =  0.0002  m,  and  the  Coriolis  parameter  /  =  1 0  4 
s^1.  The  wavefield  is  built,  as  discussed  in  Section  2.3, 
based  on  a  wind  speed  of  15  m  s_1  and  the  phase  speed 
of  the  peak  in  the  spectrum  Cp  ~  18  m  s_1.  Thus  the 
suite  of  simulations  allows  us  to  examine  a  wide  vari¬ 
ation  of  wave  age  from  swell  dominated  to  near  wind- 
wave  equilibrium.  Table  1  lists  bulk  properties  of  the 
simulations,  viz-,  the  geostrophic  wind,  wave  age,  and 
friction  velocity  u,.  U\q  is  the  reference  wind  speed  at 
a  height  of  10  m.  The  simulations  are  run  for  more  than 
50,000  timesteps  using  restart  volumes  with  fully  devel¬ 
oped  turbulence.  The  iteration  count  in  the  pressure  Pois¬ 
son  solver  is  typically  set  to  30  and  the  calculations  run 
on  either  512  or  1024  computational  cores. 


Table  1 :  Simulation  properties 


Run 

Ug  (ms  1 ) 

Cp/U  10 

m*  (m  s  *) 

A 

5 

4.8 

0.124 

B 

7.5 

3.4 

0.187 

C 

10 

2.8 

0.228 

D 

15 

1.9 

0.338 

E 

20 

1.5 

0.452 

4.  RESULTS 

Previous  field  observations  (Grachev  and  Fairall, 
2001;  Smedman  et  ah,  1999),  turbulence  closure  mod¬ 
eling  (Hanley  and  Belcher,  2008;  Makin,  2008),  and  our 
own  idealized  LES  (Sullivan  et  ah,  2008)  all  show  that 
fast  moving  swell  can  induce  marked  changes  in  the 
atmospheric  surface  layer  winds,  viz.,  an  upward  mo¬ 
mentum  flux  from  the  ocean  to  the  atmosphere,  a  low- 
level  wind  maximum,  and  departures  from  law-of-the- 
wall  scaling.  The  preliminary  LES  computations  per¬ 
formed  here  over  a  more  realistic  sea  surface  are  in  good 
qualitative  agreement  with  the  previous  studies  but  sug¬ 
gest  the  impact  of  swell  on  the  surface  layer  winds  is 
sensitive  to  the  content  of  the  wave  spectrum. 
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Figure  2:  A  snapshot  of  the  wavefield  height  h(x,y,t ) 
that  is  imposed  at  the  bottom  of  the  LES  code,  h  is 
built  from  a  sum  of  linear  plane  waves  (7).  Waves  prop¬ 
agate  left  to  right  according  to  the  dispersion  relation¬ 
ship.  The  horizontal  grid  spacing  matches  the  LES,  i.e., 
Ax  =  Ay  =  2.5  m.  The  color  bar  is  in  units  of  meters. 


One  of  the  surprising  results  from  the  present  simula¬ 
tions  is  the  significant  impact  of  swell  on  the  coherence 
and  magnitude  of  the  near-surface  pressure  fluctuations. 
This  is  illustrated  in  figure  4  where  we  compare  p' / p  for 
two  levels  of  wind  forcing  f/?  =  (5,20)  ms"1,  i.e.,  a  low- 
wind  situation  with  swell  and  a  high  wind  case  approach¬ 
ing  wind-wave  equilibrium.  The  difference  in  the  pres¬ 
sure  signals  is  striking  and  even  more  remarkable  in  ani¬ 
mations  of  the  pressure  field.  In  the  low-wind  swell  case 
there  is  a  very  strong  correlation  between  p' / p  <  0  and 
wave  crests  and  similarly  between  p' / p  >  0  and  wave 
troughs  that  extends  over  the  depth  of  the  surface  layer. 
Inspection  of  the  flow  visualization  and  animations  re¬ 
veals  that  the  strong  correlation  persists  across  the  range 
of  resolved  waves,  i.e.,  both  large  and  small  scale  waves 
appear  to  induce  a  similar  pressure  pattern.  The  coher¬ 
ence  of  the  wave  induced  pressure  field  can  extend  to  20 
m  or  more  depending  on  the  amplitude  of  the  underly¬ 
ing  wave.  Also,  the  pressure  signatures  propagate  at  the 
speed  of  the  wavefield,  additional  evidence  that  the  sig¬ 
nals  are  generated  by  surface  waves  and  not  atmospheric 
processes.  These  are  clear  signatures  of  “wave  pumping” 
by  the  surface  wavefield  on  the  atmosphere.  The  ampli- 
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Figure  3:  An  instantaneous  x  —  z  slice  of  the  3D  time  varying  computational  mesh  in  the  lowest  portion  of  the  PBL. 
The  (4,r|)  gridlines  become  level  surfaces  at  about  100  m  above  the  water.  Only  a  fraction  of  the  grid  is  displayed. 


tude  of  the  wave  spectrum  (and  hence  the  level  of  wave 
forcing)  is  held  constant  in  our  simulations  but  the  mag¬ 
nitude  of  the  turbulence,  as  measured  by  ut,  increases 
substantially  with  increasing  wind  speed.  The  structure 
of  the  near  surface  pressure  field  is  a  result  of  these  two 
competing  effects.  At  low  winds  coherent  pressure  sig¬ 
nals  are  generated  by  the  wave  motions  when  the  turbu¬ 
lence  is  weak  but  this  coherence  is  destroyed  by  strong 
turbulence  at  higher  winds. 

Figure  5  shows  that  the  impact  of  wave  age  also  ap¬ 
pears  in  the  vertical  velocity  fields.  In  the  low-wind  swell 
regime  we  observe  large-amplitude  large-scale  fluctua¬ 
tions  in  w' .  At  higher  winds  the  spatial  coherence  of  W 
is  destroyed  by  strong  turbulence.  Note  each  panel  in 
figure  5  is  sampled  at  the  sample  height  above  the  wave- 
field.  Also,  the  fields  are  made  dimensionless  by  friction 
velocity  w*  which  further  illustrates  the  strong  impact  of 
the  wave  motions  on  the  winds  in  the  surface  layer. 

In  figure  6  we  compare  vertical  profiles  of  the  mean 
wind  speed  and  turbulence  variances  for  the  different 
simulations.  These  statistics  are  computed  by  averag¬ 
ing  in  computational  coordinates,  i.e.,  across  horizontal 
planes  at  constant  vertical  height  C  Similar  to  Sullivan 
et  al.  (2008)  we  find  that  the  wind  speed  and  turbulence 
variances  depend  on  wave  age.  At  high  winds  as  the 
simulations  approach  wind-wave  equilibrium,  the  non- 
dimensional  wind  profile  (U) /u,  smoothly  approaches 
the  variation  predicted  by  law-of-the-wall.  Significant 
differences  are  observed  for  the  cases  dominated  by 
swell:  the  surface  layer  winds  are  accelerated  compared 
to  rough  wall  scaling.  As  suggested  by  the  flow  visu¬ 
alization,  the  turbulence  variances  respond  to  the  wave 
motion  in  dramatic  ways.  The  horizontal  and  vertical 
variances  are  significantly  enhanced  by  the  motion  of  the 


wave  surface  in  the  low-wind  cases.  Even  though  the  tur¬ 
bulence  is  relatively  weak  the  turbulence  variances  are 
large  near  the  wave  surface  due  to  wave  pumping. 

5.  SUMMARY 

A  large-eddy  simulation  (LES)  model  for  the  marine 
atmospheric  planetary  boundary  layer  (PBL)  is  coupled 
to  a  3D  time-dependent  surface  gravity  wavefield.  A  co¬ 
ordinate  transform  from  physical  to  computational  space 
is  used  that  accounts  for  vertical  movement  of  the  mesh 
in  physical  space.  We  use  the  geometric  conservation 
law  (GCL)  (Thomas  and  Lombard,  1979)  to  solve  for  the 
grid  speeds  that  enter  into  the  advection  of  momentum 
and  scalars.  The  algorithm  is  used  to  carry  out  a  se¬ 
ries  of  simulations  over  a  broadband  moving  wavefield 
that  conforms  to  a  Pierson-Moskowitz  wave  spectrum. 
The  wave  age  Cp/U\q  =  [1 .5,4.8]  varies  from  near  wind- 
wave  equilibrium  to  a  low-wind  swell  dominated  regime. 
In  the  low  wind  cases  we  find  features  similar  to  previous 
observational  and  modeling  investigations:  the  surface 
layer  winds  show  clear  departures  from  rough  wall  law- 
of-the-wall  scaling.  The  coherence  and  magnitude  of  the 
pressure  field  p' / p  depends  critically  on  the  motion  of 
the  underlying  wavefield  and  the  turbulence  level. 
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Figure  4:  Snapshot  of  static  pressure  fluctuations  p'  / p  in  an  x  —  z  plane  near  the  water  surface.  The  upper  panel  is  a 
swell  dominated  regime  with  wave  age  ~  4.8  while  the  lower  panel  is  a  case  near  wind-wave  equilibrium  with  wave 
age  ~  1 .4.  The  wave  spectrum  is  a  Pierson-Moskowitz  spectrum.  Notice  the  coherence  between  the  wave  field  and 
the  pressure  fluctuations  in  the  case  with  swell.  The  color  bar  is  in  units  of  m  s'  2  and  the  range  is  different  between 
the  two  cases. 
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Figure  5:  Snapshot  of  resolved  vertical  velocity  fluctuations  w'/m*  in  a  wave  following  x—y  plane  near  the  water 
surface  £  =  2.5  m.  The  left  panel  is  a  swell  dominated  regime  with  wave  age  ~  4.8  while  the  right  panel  is  a  case  near 
wind-wave  equilibrium  with  wave  age  ~  1 .4.  The  wave  spectrum  at  the  bottom  of  the  PBLs  is  the  same.  Notice  the 
range  of  the  color  bar  is  different  between  the  two  cases.  The  (normalized)  fluctuations  in  the  wind-wave  equilibrium 
case  are  smaller. 
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Figure  6:  Vertical  profiles  of  wind  speed  (left  panel)  and  turbulence  variances  (right  panels)  for  different  values 
of  wave  age  Cp/U\q.  Friction  velocity  n*  is  used  for  normalization.  The  dashed  black  line  is  the  rough  wall  formula 
U /m*  =  In {z/zo) /k,  where  K  =  0.4.  Temporal  and  spatial  averaging  is  used  to  make  the  statistics.  The  spatial  averaging 
is  over  ^  —  1]  planes,  i.e.,  at  constant  C 
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1.  INTRODUCTION 

A  wide  spectrum  of  atmospheric  motions  impacts  the 
flow  environment  in  and  around  wind  farms.  The  im¬ 
portant  time  and  space  scales  span  an  enormous  range, 
large-scale  long-time  decadal  climate  down  to  rapidly 
evolving  thin  boundary  layers  that  grow  on  an  individual 
turbine  blade  (Schreck  et  al.,  2008).  In  this  scale  hier¬ 
archy,  the  atmospheric  planetary  boundary  layer  (PBL), 
where  the  important  scales  range  from  O(m)  or  less  to 
O(km)  or  more,  plays  a  critical  if  not  dominant  role  in 
setting  the  performance  of  an  individual  turbine  as  well 
as  the  power  output  of  an  entire  wind  park.  PBL  mo¬ 
tions  are  turbulent  (three-dimensional  and  time  depen¬ 
dent)  and  couple  to  larger-scale  atmospheric  motions  and 
land  use,  e.g.,  weather  events  diurnal  heating  and  cool¬ 
ing,  thermal  stratification,  surface  roughness,  vegetative 
canopies,  wind  waves  and  local  orography  all  influence 
wind  turbine  performance  to  varying  degrees.  For  exam¬ 
ple,  the  afternoon  collapse  of  the  heated  daytime  PBL 
over  the  Great  Plains  followed  by  surface  cooling  can 
lead  to  a  windy  weakly  stratified  boundary  layer  with 
a  nocturnal  jet  positioned  near  hub  height  z  ~  100  m 
(Beare  et  al.,  2006;  Banta  et  al.,  2008).  Turbulence  in 
this  regime  is  episodic,  non-Gaussian,  and  can  interact 
with  gravity  waves  and  trigger  Kelvin-Helmholtz  insta¬ 
bilities  leading  to  intermittent  loads  that  fatigue  turbine 
components  (Kelley  et  al.,  2003). 

The  objective  of  the  present  work  is  to  describe  our  re¬ 
cent  developments  in  constructing  and  utilizing  a  large- 
eddy  simulation  (LES)  model  for  the  atmospheric  PBL 
where  the  lower  boundary  shape  is  modestly  complex, 
i.e.,  we  define  modestly  complex  as  orography  of  height 
h  to  be  a  single  valued  function  of  the  horizontal  coor¬ 
dinates  h  =  h(x.y).  This  enhanced  simulation  capabil¬ 
ity  will  allow  us  to  examine  basic  interactions  between 
stratified  PBL  turbulence  and  landscape  features,  but  also 
provide  information  about  the  winds  turbines  might  be 
exposed  to  in  a  local  undulating  environment.  Under¬ 
standing  the  flow  environment  created  by  complex  small- 
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Figure  1:  Cathedral  Rocks  wind  farm  located 

above  an  escarpment  on  the  Eyre  Peninsula  South 
Australia  adapted  from  http://www.yorkcivil.com.au/ 
projects/projects/49/cathedral_rocks_wind_farm. 


scale  topography,  as  in  figure  1,  is  of  particular  impor¬ 
tance  because  steep  slopes  often  generate  high  levels  of 
turbulence  in  zones  of  intermittent  flow  separation.  Cou¬ 
pled  with  background  PBL  turbulence  this  orographi- 
cally  generated  turbulence  can  significantly  restrict  the 
area  available  for  turbine  placement,  in  which  the  flow 
conforms  to  turbine  design  standards  (IEC,  2005;  Ayotte, 
2008).  We  emphasize  that  a  complete  simulation  of  all 
the  turbulence  length  and  time  scales  generated  by  the 
terrain  in  figure  1  far  exceeds  current  computational  ca¬ 
pabilities. 

2.  LES  ALGORITHM  WITH  SURFACE  TER¬ 
RAIN 

2.1  LES  with  a  flat  lower  boundary 

Typical  LES  model  equations  for  a  dry  Boussinesq 
PBL  include  at  a  minimum:  a)  transport  equations  for 
momentum  pu;  b)  a  transport  equation  for  a  conserved 
buoyancy  variable  (e.g.,  virtual  potential  temperature  0); 
c)  a  discrete  Poisson  equation  for  a  pressure  variable 
II  to  enforce  incompressibility;  and  closure  expressions 
for  subgrid-scale  (SGS)  variables,  e.g.,  a  subgrid-scale 


equation  for  turbulent  kinetic  energy  e.  Our  notation 
denotes  a  resolved  scale  variable  by  an  overbar  ( )  and 
thus  u  are  spatially  filtered  Cartesian  velocity  compo¬ 
nents.  In  our  LES  code  with  a  flat  boundary  the  spatial 
discretization  is  second-order  finite  difference  in  the  ver¬ 
tical  z  direction  and  pseudospectral  in  horizontal  x  —  y 
planes.  Thus  a  staggered  arrangement  of  variables  is 
used  with  (17,  v.  n.  0)  stored  at  cell  centers  and  (w.e)  lo¬ 
cated  at  cell  faces;  this  variable  layout  is  advantageous 
because  it  tightly  couples  velocity  and  pressure  in  our  in¬ 
compressible  formulation.  The  equations  are  integrated 
forward  in  time  using  a  fractional  step  method  utilizing 
dynamic  time  stepping  with  a  fixed  Courant-Fredrichs- 
Lewy  (CFL)  number  (Sullivan  et  al.,  1996;  Spalart  et  al., 
1991).  The  code  parallelization  is  accomplished  us¬ 
ing  the  Message  Passing  Interface  (MPI)  and  a  2D  do¬ 
main  decomposition.  Simulations  have  used  as  many  as 
16,384  computational  cores  for  meshes  with  30723  grid- 
points  (Sullivan  and  Patton,  2010). 


2.2  Coordinate  transformation 

In  order  to  adapt  the  LES  model  with  a  flat  lower  bot¬ 
tom,  outlined  in  Section  2.1,  to  an  atmospheric  PBL  flow 
with  a  varying  boundary  shape  we  apply  a  conventional 
grid  transformation  to  the  LES  equations.  The  transfor¬ 
mation  maps  the  surface  following  non-orthogonal  mesh 
onto  a  flat  computational  space  according  to  the 

rule: 
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The  Jacobian,  which  is  needed  to  move  between  physical 
and  computational  spaces,  is  (Anderson  et  al.,  1984,  see 
Eq.  5-234) 
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Several  strategies  are  available  for  building  the  com¬ 
putational  mesh  in  physical  space.  Most  often  we  em¬ 
ploy  simple  algebraic  stretching  (e.g.,  Anderson  et  al., 
1984,  see  page  358).  This  technique  builds  a  smoothly 
varying  mesh  along  each  vertical  coordinate  line.  The 
stretching  factor  K  =  is  the  grid  spacing  ra¬ 

tio  between  neighboring  k  and  k  +  1  vertical  gridpoints. 
A  typical  value  is  K  ~  1.03  or  less  when  A Z\  =  1  m  and 
the  top  of  the  domain  is  1000  m  and  the  number  of  verti¬ 
cal  gridpoints  N-  =  256. 


2.3  LES  equations  in  curvilinear  coordinates 

The  LES  equations  in  curvilinear  coordinates  can  be 
derived  in  a  straightforward  fashion  by  applying  the 
chain  rule  for  differentiation.  The  set  of  LES  equations 
written  in  computational  coordinates  under  the  transfor¬ 
mation  (1)  and  (2)  are: 
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The  equations  are  expressed  in  strong  conservation  form 
using  the  “contravariant  flux”  velocity 
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(3a)  is  the  mass  conservation  (continuity)  equation,  (3b) 
is  the  momentum  transport  equation,  (3c)  is  the  scalar 
transport  equation,  (3d)  is  the  subgrid-scale  energy  trans¬ 
port  equation,  and  (3e)  is  the  pressure  Poisson  equation. 
The  right  hand  sides  of  (3)  model  physical  processes  in 
the  atmospheric  PBL,  e.g.,  pressure  gradients,  Coriolis 
rotation,  divergence  of  subgrid-scale  fluxes,  buoyancy, 
and  in  the  case  of  the  SGS  e  equation  also  diffusion  and 
dissipation. 

A  key  step  in  this  formulation  is  co-locating  the  solu¬ 
tion  variables  (u,n,0,e)  at  cell  centers  (Sullivan  et  al., 
2000,  2008)  which  leads  to  a  simple  compact  differenc¬ 
ing  stencil.  In  order  to  maintain  tight  velocity-pressure 
coupling  the  flux  velocities  (U  .  V)  are  located  at  cell 
centers  while  W  is  located  at  a  cell  face.  This  mim¬ 
ics  the  variable  layout  in  our  usual  staggered  code  with 
a  flat  bottom.  As  in  Sullivan  et  al.  (2008)  we  use  mo¬ 
mentum  interpolation  to  construct  the  intermediate  right 
hand  sides  for  the  flux  velocities.  Formally,  the  equa¬ 
tion  set  (3)  has  the  same  structure  as  in  the  case  with  a 
flat  lower  boundary  and  thus  similar  spatial  and  temporal 
discretizations  are  used  to  advance  them  forward  in  time. 
Thus  the  spatial  differencing  is  pseudospectral  in  the  hor¬ 
izontal  computational  directions  fq.  q  )  and  second-order 
finite  differences  in  the  ^-coordinate. 

The  major  algorithmic  change,  compared  to  the  flat 
LES  code,  is  the  pressure  formulation.  As  a  conse¬ 
quence  of  the  incompressible  flow  assumption  and  the 
non-orthogonal  mesh  we  are  forced  to  solve  a  general 


Poisson  equation  (3e)  for  pressure.  A  direct  solution  of 
(3e)  is  not  possible.  Inspection  of  (3a)  and  (3e)  suggests 
the  stationary  iteration  scheme 

®(n'+1)  =  z>( nf)  -  (5) 

for  II  where  At  is  the  timestep  the  superscript  i  is  an  iter¬ 
ation  index  and  the  operator  '!)  is  an  approximate  diago- 
nalization  of  the  operator  At:  At  is  the  complete  left  hand 
side  of  (3e).  At  convergence  £>( IT)  =  £>(IT+1)  and  (5) 
numerically  satisfies  mass  conservation.  (5)  is  designed 
so  that  it  can  be  easily  inverted  using  a  combination  of 
2D  Fast  Fourier  Transforms  and  tridiagonal  matrix  in¬ 
versions.  Furthermore  it  nicely  maps  onto  our  2D  MPI 
parallelization. 

The  algorithm  outlined  above  is  a  significant  advance 
over  our  previous  scheme  (Sullivan  et  al.,  2000,  2008). 
It  uses  a  more  general  grid  transformation  that  permits 
3D  lower  boundary  shapes,  it  has  a  more  general  and  ro¬ 
bust  pressure  solver,  and  it  has  an  improved  treatment  of 
the  surface  boundary  conditions.  We  have  used  as  many 
as  4096  computational  cores  for  high  resolution  simu¬ 
lations,  (1024  x  384  x  256)  gridpoints,  with  the  above 
scheme.  Further  details  of  the  algorithm  including  the 
posing  of  the  rough  wall  boundary  conditions  will  be  de¬ 
scribed  in  a  future  publication. 

3.  SIMULATION  STRATEGY 

A  series  of  LES  are  performed  to  highlight  the  in¬ 
teractions  between  small-scale  terrain  and  PBL  turbu¬ 
lence  and  to  exercise  the  code  for  both  two-dimensional 
and  three-dimensional  surface  shapes.  In  these  compu¬ 
tations,  we  simplify  the  external  forcing  compared  to  an 
atmospheric  PBL,  i.e.,  we  drive  the  flows  with  a  con¬ 
stant  large-scale  pressure  gradient  Tx/ p  oriented  along 
the  x-  direction  and  turn  off  buoyancy  influences.  Hence 
the  simulations  are  similar  in  spirit  to  a  wind  tunnel  con¬ 
figuration.  The  lower  surface  is  assumed  to  be  fully 
rough  and  we  impose  z0  boundary  conditions  based  on 
the  winds  at  the  first  gridpoint  off  the  surface  ( e.g Mo- 
eng,  1984). 

The  simulations  with  terrain  are  computationally  ex¬ 
pensive  compared  to  flat  wall  LES  especially  when  the 
slopes  of  the  boundary  shape  are  steep.  In  these  simu¬ 
lations,  20  to  30  iterations  of  (3e)  are  required  to  reduce 
the  residual  of  the  continuity  equation  (3a)  to  near  ma¬ 
chine  accuracy.  This  increases  the  computational  time  by 
a  factor  of  two  or  more  compared  to  a  flat  wall  simula¬ 
tion.  In  addition,  we  use  fine  mesh  resolution  O(m)  (see 
below)  and  thus  the  timesteps  are  pushed  to  small  values 
by  the  CFL  constraint  in  our  shear  driven  computations. 
Le  and  Moin  (1991)  and  Zheng  and  Petzold  (2006)  pro¬ 
pose  fractional  step  schemes  that  require  only  a  single 


pressure  projection  step  in  their  multi-stage  Runge-Kutta 
time  stepping  schemes  which  would  potentially  lower  the 
cost  of  including  terrain  in  an  incompressible  LES. 

To  allow  efficient  simulations  of  the  various  flows  de¬ 
scribed  in  Section  4.2  we  adopt  the  following  strategy: 
1)  the  simulations  are  first  carried  forward  with  a  flat  bot¬ 
tom  until  the  turbulence  is  fully  developed  and  the  hor¬ 
izontally  averaged  momentum  flux  profile  is  linear  in  z 
over  the  computational  domain;  2)  we  then  restart  the 
simulation  with  the  curved  lower  bottom  using  a  data 
volume  archived  at  late  time  from  the  flat  wall  case. 
These  simulations  are  then  advanced  in  time  until  the 
turbulence  is  nearly  re-cycled  through  the  computational 
box.  Thus  our  strategy  mimics  the  conditions  in  a  wind 
tunnel  where  a  boundary-layer  first  develops  over  a  flat 
wall  and  then  encounters  an  obstacle  far  downstream 
of  the  inlet.  This  simulation  technique  is  hinted  at  by 
Gong  et  al.  (1996,  see  discussion  starting  on  page  24) 
and  allows  both  spatial  averaging  over  homogeneous  di¬ 
rections  (e.g.,  the  y— direction)  and  over  a  set  of  realiza¬ 
tions.  Multiple  realizations  can  be  created  by  simply  re¬ 
starting  the  terrain  simulations  using  different  initial  vol¬ 
umes  from  the  flat  wall  case.  The  transient  induced  by  in¬ 
serting  the  bottom  terrain  is  short  lived  and  not  analyzed. 
The  pressure  solver  is  able  to  generate  incompressible 
flow  in  a  single  timestep  after  inserting  the  terrain  since 
it  is  developed  from  the  discrete  version  of  the  continuity 
equation. 

4.  SAMPLE  RESULTS 

4.1  Turbulent  flow  past  two-dimensional  shapes 

Taylor  (1998)  provides  the  variation  of  form  (pressure) 
drag  versus  waveslope  ak  for  linearized  mean  flow  mod¬ 
els  with  different  turbulent  closures.  The  bumps  are  two- 
dimensional  (no  y  variation)  with  shape  h  =  a  cos  2kx/X: 
the  wavelength  X  and  wavenumber  k  are  related  by  X  = 
2%/k.  We  re -plot  their  results  in  figure  2.  At  low  waves¬ 
lope,  ak  <  0.17,  the  theoretical  results  are  nearly  inde¬ 
pendent  of  the  closure.  Above  ak  >  0.3  the  bumps  are 
steep  and  the  linearized  calculations  tend  to  breakdown. 

We  perform  3D  LES  of  a  similar  turbulent  boundary- 
layer  flow  over  2D  sinusoidal  bumps  using  a  modest 
grid  mesh  of  (256,256, 128)  gridpoints  with  a  relatively 
small  surface  roughness  X/z0  —  5  x  105.  The  computa¬ 
tions  are  carried  out  for  four  different  waveslopes  ak  = 
(0.1,0.25,0.35,0.5).  At  low  ak  the  LES  values  closely 
match  the  linearized  calculations.  The  LES  continues 
to  work  smoothly  for  waveslopes  as  large  as  0.5  (this 
is  the  largest  value  tested).  An  additional  LES  with  a 
relatively  large  surface  roughness  X/z0  =  1  x  103  is  also 
shown  in  figure  2.  This  value  of  roughness  causes  large 
flow  separation  with  flow  reattachment  on  the  forward 


face  of  the  upstream  wave.  This  is  clearly  observed  in 
the  visualization  of  the  pressure  contours  and  flow  vec¬ 
tors.  It  is  interesting  that  the  pressure  drag  for  the  small 
and  large  roughness  are  almost  identical.  We  find  that  in 
the  large  roughness  case  the  flow  above  the  bumps  tends 
to  skip  from  crest-to-crest  with  relatively  slow  recircu¬ 
lating  flow  in  the  wave  troughs.  In  the  small  roughness 
(non-separated)  case  the  wave  signature  is  clearly  visible 
in  the  vertical  velocity  field  but  is  destroyed  by  vigor¬ 
ous  turbulence  generated  by  flow  separation  in  the  large 
roughness  case  (see  figure  2). 

4.2  Turbulent  flow  past  three-dimensional  shapes 

Hills,  ridges,  bluffs,  land-sea  escarpments,  etc.  tend  to 
generate  a  local  speedup  in  the  boundary-layer  winds  and 
thus  are  potential  targets  for  wind  turbine  sites  (see  fig¬ 
ure  1).  Often  these  orographic  features  are  geometrically 
complex  and  3D,  i.e.,  with  their  characteristic  horizontal 
lengths  and  widths  of  similar  scale.  Thus  it  is  important 
to  examine  the  structure  of  the  boundary-layer  winds  in 
the  presence  of  3D  obstacles. 

We  compute  turbulent  flow  over  and  around  three 
canonical  shapes,  viz.,  a  hill,  gap,  and  crater.  This 
demonstrates  the  LES  code’s  ability  to  handle  mod¬ 
estly  complex  3D  orography  and  further  illustrates  the 
rich  level  of  fluid  dynamical  phenomenon  generated  by 
the  interactions  between  boundary-layer  turbulence  and 
small-scale  landscape  features.  In  the  following  sec¬ 
tions,  the  boundary-layer  flows  are  created  using  the 
strategy  described  in  Section  3.  The  simulation  with  a 
flat  wall  is  first  run  for  about  130,000  timesteps  which 
is  ~  7  large  eddy  turnover  times.  Each  simulation  uses 
a  unique  lower-bottom  shape  and  generates  a  terrain  fol¬ 
lowing  grid  using  smooth  vertical  grid  stretching.  The 
computational  box  (. Lx,Ly,Lz )  =  (2560,640, 1000)  m  the 
number  of  gridpoints  in  the  three  coordinate  directions 
(Nx.Ny,Nz)  =  (1024,256,256)  and  the  horizontal  spac¬ 
ing  (Ax,  Ay)  =  (2.5, 2.5)  m.  The  vertical  spacing  is  1 
m  at  the  surface  and  smoothly  increases  to  the  top  of  the 
box.  Approximately  40  gridpoints  are  located  in  the  first 
50  m  above  the  surface.  The  large-scale  pressure  gra¬ 
dient  is  set  to  !?y/p  =  1.63  10  4  m  s“2  which  generates 
winds  of  about  10.5  m  s~*  at  the  top  of  the  box.  The  sur¬ 
face  roughness  z0  =  0.05  m.  The  flow  blockage  is  small 
since  h/Lz  <  0.05  everywhere. 

4.2.1  Isolated  hill 

An  isolated  3D  hill  (see  figure  3)  is  positioned  in 
the  simulation  with  its  summit  located  at  (xc,yc)  = 
(1280,320)  m.  Geometrical  properties  of  the  cosine 
shaped  hill  are:  summit  height  of  50  m,  maximum  slope 
of  0.6,  characteristic  length  scale  L  =  67  m,  and  the  max¬ 


imum  x  —  y  planform  is  approximately  a  circle  of  diame¬ 
ter  equal  to  4A  =  268  m.  The  shape  of  the  hill  is 

h(Ay')  =  t(1+cos^;)  (1+cos^r)  (6) 

where  b  =  25  m,  x'  =  x  —  xc,  and  (|xy|,  |y'|)  <  2L. 

Figure  3  shows  an  instantaneous  snapshot  of  typical 
flow  features  that  develop  around  the  steep  isolated  hill. 
The  overall  impression  is  that  these  flow  patterns  are 
unique  compared  to  a  2D  case  with  similar  characteris¬ 
tic  length  scale  L  and  maximum  slope.  First,  a  complex 
flow  separation  pattern  is  generated  -  there  is  flow  sep¬ 
aration  over  the  hill  crest  as  in  the  2D  case  but  separa¬ 
tion  also  occurs  along  the  flanks  of  the  hill.  The  region 
of  intense  low  pressure  along  the  hill  crest  is  confined 
to  a  distance  of  about  2. 5 A  in  the  y— direction.  Anima¬ 
tions  show  that  the  spanwise  pressure  pattern  along  the 
line  x 1  —  0  is  occasionally  interrupted  by  tongues  of  low 
pressure  that  arc  around  the  hill  sides  in  a  crescent  shape. 

The  most  striking  feature  of  this  simulation  however 
is  the  complex  multi-scale  wake  flow  that  develops  aft 
of  the  hill.  It  is  a  collection  of  temporally  evolving  vor¬ 
tices  with  their  primary  axis  aligned  with  z.  At  any  time 
there  can  be  multiple  vortices  of  various  scales  present 
but  often  two  larger  vortices  dominate  the  wake  as  shown 
in  figure  3.  There  is  a  strong  return  flow  along  the  line 
y'  =  0  separating  the  two  vortices.  Notice  also  that  the 
location  of  the  vortices  is  well  downstream  of  the  hill 
summit  near  the  boundary  where  the  hill  flattens  out  and 
blends  into  the  flat  surface;  this  is  clearly  downwind  of 
the  region  of  maximum  hill  slope.  Note  we  are  using  low 
pressure  as  a  criterion  for  vortex  identification,  see  Lin 
et  al.  (1996)  for  a  comparison  of  methods.  The  coher¬ 
ent  rotation  of  the  surface  velocity  vectors  is  additional 
evidence  that  the  regions  of  low  pressure  are  indeed  vor¬ 
tical  cores.  Broadly,  the  mean  flow  patterns  shown  in 
figure  3  are  similar  to  the  early  stage  of  a  flow  around  a 
barchan  sand  dune  described  by  Ortiz  and  Smolarkiewicz 
(2009),  and  the  flow  measurements  around  an  axisym- 
metric  bump  reported  by  Byun  and  Simpson  (2006). 


4.2.2  Gap  flow 

Flows  in  gaps  separating  ridges  and  hills  are  a  com¬ 
mon  landscape  feature.  The  ability  of  the  LES  to  simu¬ 
late  this  type  of  turbulent  flow  is  shown  in  figure  4.  The 
shape  of  the  terrain  is 


(7) 


where  the  spanwise  extent  and  depth  of  the  gap  are  con¬ 
trolled  by  the  function 
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Figure  4  illustrates  the  speedup  of  the  u—  compo¬ 
nent  of  the  horizontal  wind  in  a  narrow  gap.  Notice  the 
speedup  is  largest  in  the  region  where  the  ridge  begins  to 
blend  into  the  gap  floor,  i.e.,  along  the  sides  of  the  gap. 
Visualization  of  the  pressure  field  shows  a  similar  feature 
where  the  low  pressure  contours  are  most  negative  on  the 
slope  approaching  the  valley  and  then  become  less  neg¬ 
ative  away  from  the  gap.  The  pattern  is  asymmetric  in 
time  shifting  from  side-to-side  in  the  gap.  This  is  slightly 
surprising  since  the  initial  expectation  is  that  the  winds 
along  the  gap  centerline  would  be  highest.  The  vertical 
y  —  z  planes  illustrate  the  vigorous  flow  separation  aft  of 
the  two  ridges. 

4.2.3  Crater  flow 

Boundary-layer  flows  in  and  above  open  and  closed 
basins  are  of  importance  for  a  variety  of  applications, 
e.g.,  diffusion  of  pollutants.  The  recent  METCRAX  field 
campaign  (Whiteman  and  CoAuthors,  2008)  focused  on 
the  complex  flow  patterns  generated  by  stably-stratified 
flow  in  a  closed  basin.  In  figure  5  we  show  an  example  of 
the  LES  code’s  capability  to  simulate  neutrally  stratified 
flow  in  an  idealized  closed  crater  whose  shape  h  is  given 
by  the  negative  of  (6).  The  nominal  diameter  is  4 L  =  400 
m,  the  depth  equals  50  m,  and  the  maximum  slope  equals 
0.39.  The  crater  center  (. xc,yc )  =  (1280, 320)  m. 

The  instantaneous  snapshot  of  the  flowfield  in  fig¬ 
ure  5  shows  several  interesting  and  complex  features. 
The  negative  (low)  pressure  contours  are  indicators  of 
small-scale  vortices  located  near  the  crater.  Animations 
show  a  rapid  evolution  of  these  vortical  structures  in 
the  crater  interior.  The  high  pressure  contours  are  in¬ 
dicators  of  the  fluctuating  re-attachment  along  the  crater 
backwall  and  the  very  negative  pressure  contours  around 
the  crater  rim  are  due  to  flow  accelerations  around  the 
crater  lip.  Visualization  in  x  —  z  planes  along  the  crater 
centerline  show  intermittent  ejections  of  fluid  into  the 
overlying  boundary-layer  flow.  We  mention  that  strong 
convergence  of  surface  streamlines,  e.g.,  as  observed  at 
(x,y)  =  1350,250)  m  in  figure  5,  is  often  an  indicator  of 
3D  flow  separation  (e.g.,  Byun  and  Simpson,  2006).  The 
flow  patterns  inside  the  crater  are  clearly  distinct  from 
those  generated  behind  an  isolated  hill  in  figure  3. 

5.  SUMMARY 

A  massively  parallel  algorithm  and  code  for  large- 
eddy  simulation  (LES)  of  atmospheric  planetary  bound¬ 


ary  layers  (PBLs)  with  modestly  complex  orography  is 
described.  Our  LES  model  equations  adopt  an  incom¬ 
pressible  Boussinesq  flow  model  with  high  Reynolds 
number  rough  wall  boundary  conditions  along  the  lower 
boundary.  A  co-located  variable  layout  and  a  conven¬ 
tional  coordinate  transformation  from  physical  to  com¬ 
putational  space  are  used.  The  grid  mesh  in  physical 
space  is  terrain  following  and  non-orthogonal,  a  more 
general  formulation  can  be  incorporated  into  the  scheme. 
The  key  new  step  compared  to  a  flat  wall  code  is  the 
formulation  of  the  pressure  equation  and  designing  an 
algorithm  for  the  solution  of  the  pressure  Poisson  equa¬ 
tion.  The  algorithm  is  sufficiently  general  to  allow  simu¬ 
lations  of  PBLs  over  a  spectrum  of  time  dependent  mov¬ 
ing  water  waves  (Sullivan  et  ah,  2010).  We  present  sev¬ 
eral  sample  calculations  of  neutrally-stratified  turbulent 
flow,  (similar  to  a  wind-tunnel  setup)  past  2D  sinusoidal 
bumps  and  3D  obstacles,  viz.,  a  hill,  gap,  and  crater. 
These  calculations  highlight  the  importance  of  flow  sep¬ 
aration  and  coupling  with  background  PBL  turbulence, 
and  the  ability  of  the  algorithm  to  simulate  turbulent 
flows  with  an  undulating  lower  boundary. 

In  the  future  we  plan  to  implement  an  algebraic 
stress  closure  model  for  subgrid-scale  fluxes  and  vari¬ 
ances  (Wyngaard,  2004)  and  validate  the  code  against 
the  wind-tunnel  measurements  of  Ayotte  and  Hughes 
(2004),  and  Gong  et  al.  (1996),  and  field  observations. 
Also,  PBL  simulations  will  be  carried  out  with  unsta¬ 
ble  and  stable  stratification.  The  fine  mesh  large-eddy 
simulations  described  here  can  be  used  in  two  ways:  1) 
they  provide  insight  about  the  fundamental  interactions 
between  turbulence  and  terrain  which  can  impact  iso¬ 
lated  wind  turbines  and  wind  parks;  and  2)  the  detailed 
datasets  can  be  used  for  building  parameterizations  of 
separated  flows. 
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Figure  2:  Turbulent  flow  over  2D  sinusoidal  bumps.  The  upper  panel  shows  the  pressure  drag  coefficient  Co  as 
function  of  waveslope.  Co  is  obtained  from  integration  of  the  pressure-waveslope  correlation  and  is  normalized  by  the 
friction  velocity  squared  n*.  The  lines  are  theoretical  predictions  from  linearized  calculations  with  different  turbulence 
closures  at  a  small  roughness  X/z0  =  104  from  Taylor  (1998).  The  orange  solid  bullets  are  results  from  LES  with  a 
similar  small  surface  roughness  X/z0  =  5  x  10s.  The  blue  triangle  is  an  LES  with  large  surface  roughness  X/z0  = 
1  x  103  which  leads  to  extensive  flow  separation  between  the  wave  crests.  The  lower  panels  are  visualization  of 
vertical  velocity  w  at  a  nominal  height  of  z  =  10  m  above  the  bumps  at  waveslope  ak  =  0.5  from  LES.  The  left  and 
right  panels  are  small  and  large  roughness  X/z0  =  (5  x  105, 1  x  103),  respectively.  The  color  bar  is  in  units  of  m  s  1 
and  is  different  between  the  two  plots. 
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Figure  3:  Turbulent  flow  around  a  steep  3D  hill.  The  upper  left  panel  shows  an  oblique  view  of  the  hill  with  the 
primary  flow  direction  parallel  to  the  x— direction.  In  the  upper  right  panel  we  show  the  time  averaged  streamline 
patterns  around  the  hill  at  the  height  z  =  5.6  m.  In  the  lower  panel  we  show  color  contours  of  fluctuating  pressure  p'  / p 
overlayed  with  horizontal  flow  vectors  at  a  nominal  height  of  z  =  5.6  m  above  the  hill  surface.  The  color  bar  is  in  units 
of  m  s'  2.  The  planform  of  the  hill  (i.e.,  the  location  where  the  cosine  shaped  hill  blends  into  the  flat  bottom  boundary) 
is  approximately  indicated  by  the  circular  white  line.  The  hill  summit  h  =  50  m  is  located  at  (xc . yc )  —  (1280,320)  m. 


Figure  4:  Turbulent  flow  in  a  gap  between  two  steep  ridges.  The  primary  flow  direction  is  parallel  to  the  x—  direction. 
Contours  of  the  instantaneous  u— velocity  component  are  shown  and  the  color  bar  is  in  units  of  m  s  .  The  two  y  —  z 
planes  are  located  at  x  =  (1500,2000)  m,  and  the  ridgeline  is  located  at  x  =  1280  m.  The  horizontal  plane  is  about  4 
m  off  the  surface.  The  lower  150  m  of  the  boundary  layer  is  depicted. 
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Figure  5:  Neutrally  stratified  turbulent  flow  in  a  crater:  The  upper  left  panel  is  an  oblique  view  of  the  crater  geometry. 
The  upper  right  panel  shows  instantaneous  static  pressure  contours  p’ / p  and  horizontal  velocity  vectors  at  a  nominal 
height  of  z  =  2.5  m  above  the  surface.  The  white  circle  is  approximately  the  outline  of  the  crater  rim.  The  color  bar 
is  in  units  of  m  s“2.  The  bottom  panel  is  a  150  s  time  average  of  the  streamlines  and  pressure  field  in  an  x  —  z  plane 
along  the  crater  centerline. 
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Momentum  and  scalar  fluxes  at  the  air-sea  interface  depend  on  numerous  small  and  large  scale 
processes  that  couple  atmospheric  turbulence  and  the  underlying  surface  gravity  wave  field  (Sul¬ 
livan  and  McWilliams,  2010).  The  surface  wave  field  moves  randomly,  propagates  rapidly,  breaks 
intermittently,  and  supports  wind  stress.  Wind-wave  coupling  determines  the  sea  surface  drag, 
and  at  high  winds  provides  Oil)  control  on  hurricane  intensity.  Despite  intense  study,  the 
present  understanding  of  wind-wave  coupling  is  incomplete  and  hence  our  ability  to  predict  high 
wind  storm  dynamics,  such  as  hurricanes,  is  limited  (Emanuel,  2004;  Sanford  et  al.,  2011).  A 
fundamental  unresolved  observational  challenge  is  determining  the  correlation  between  surface 
pressure  and  wave  slope  (form  drag)  near  an  undulating  air-water  interface  with  breaking  waves. 

Our  past  modeling  using  direct  numerical  and  large-eddy  simulations  (DNS  and  LES)  over  an 
idealized  monochromatic  surface  wave  have  provided  insight  about  wind-wave  coupling  processes, 
e.g.,  critical  layer  dynamics,  wave  driven  winds,  the  transport  of  fluxes  in  non-equilibrium  wind- 
wave  conditions  and  the  impact  of  wave  state  on  the  marine  boundary  layer  (Sullivan  et  ah, 
2000;  Sullivan  and  McWilliams,  2002;  Sullivan  et  ah,  2008).  These  simulation  results  have  also 
aided  the  interpretation  of  observations.  The  computational  project  proposed  here  builds  on 
these  past  efforts  but  is  a  significant  advancement  as  it  seeks  to  examine  the  coupling  between 
stratified  atmospheric  turbulence  in  a  high  wind  marine  boundary  layer  and  a  spectrum  of  phase- 
resolved  surface  waves  using  newly  developed  large-eddy  simulation  technology.  The  high- wind 
flow  regime  of  interest  represents  a  severe  computational  challenge  as  it  couples  atmospheric 
turbulence  and  3D  moving  surface  wave  fields  over  a  wide  range  of  scales,  0(  1  —  1000)  m.  The 
very  high  spatial  and  temporal  resolution  of  the  proposed  simulations  will  provide  information 
about  the  impacts  of  the  wave  field  on  drag,  turbulence  fluxes  and  variances,  mean  wind  profiles, 
and  coherent  structures.  The  proposed  project  is  an  integral  component  of  the  ffigh  Resolution 
Air-Sea  Interaction  (Hi-Res)  observational  study  shown  in  Fig.  1.  Hi-Res  sampled  the  spatial 
and  temporal  evolution  of  the  wave  fields  and  near  surface  boundary  layer  extensively  at  wind 
speeds  varying  from  10-20  ms-1.  A  longer  term  objective  of  the  proposed  project  is  to  use 
simulation  data  to  compare  with  and  more  fully  interpret  observations  collected  in  Hi-Res. 

Our  computational  project  supports  ongoing  research  funded  by  the  Physical  Oceanography 
Program  at  the  Office  of  Naval  Research  and  the  Department  of  Energy.  Also,  the  proposed 
project  addresses  the  NCAR  Imperative  and  Frontiers  on  prediction  of  severe  weather  and  winds 
for  renewable  energy^,  and  the  NCAR  Wyoming  Supercomputer  Community  Science  Objectives 
“  ...  representation  of  detailed  processes  at  the  air- water  interface”,  and  boundary  layer 
modeling  for  studies  of  wind  energy.” 
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C.  Computational  plan 

1.  Numerical  method 

We  briefly  outline  the  important  aspects  of  our  numerical  scheme  with  moving  waves.  The 
governing  equations  which  describe  three-dimensional  time-dependent  turbulent  winds  in  a  dry 
incompressible  Boussinesq  atmospheric  boundary  layer  include:  a)  three  transport  equations  for 
momentum  pu;  b)  a  transport  equation  for  a  conserved  buoyancy  variable  {e.g.,  virtual  potential 
temperature  0);  c)  a  discrete  Poisson  equation  for  a  pressure  variable  II  to  enforce  incompressibil¬ 
ity;  and  closure  expressions  for  subgrid-scale  (SGS)  variables,  e.g.,  a  subgrid-scale  equation  for 
turbulent  kinetic  energy  e,  e.g.,  see  Sullivan  and  Patton  (2011).  The  physical  processes  included 
in  the  LES  boundary-layer  equations  include,  temporal  time  tendencies,  advection,  pressure  gra¬ 
dients,  Coriolis  rotation,  divergence  of  subgrid-scale  fluxes,  buoyancy,  resolved  turbulence,  and 
in  the  case  of  the  SGS  e  equation  also  diffusion  and  dissipation. 

In  order  to  simulate  turbulent  flow  over  a  multi-component  propagating  surface  wave  held, 
we  have  adapted  our  LES  model  with  a  flat  bottom  to  the  situation  with  a  three-dimensional 
time-dependent  boundary  shape  h  =  h(x,  y,  t )  by  applying  a  transformation  to  the  physical  space 
coordinates  (x,  y,  z )  that  maps  them  onto  computational  coordinates  (£,  77,  £).  The  computational 
mesh  in  physical  space  is  surface  following,  non-orthogonal,  and  time  varying.  Vertical  gridlines 
are  held  fixed  at  a  particular  {x,  y)  location  on  the  surface  but  the  lines  undergo  vertical  trans¬ 
lation  as  a  function  of  time  t,  i.e.,  vertical  gridlines  are  wave  following,  see  Figs.  2  and  3.  The 
LES  equations  are  expressed  in  strong  conservation  form  using  the  “contravariant  flux”  velocity 

tj  _  flj_ 

i  ~  J  dXj 


with  the  advective  term  for  variable  ^  evaluated  in  flux  form,  i.e.,  dUjfltjd^j.  In  (1),  ( J,  d^i/dxfl 
are  the  Jacobian  and  metric  connection  functions  linking  physical  and  computational  space. 

A  key  step  in  the  wavy  LES  formulation  is  co-locating  the  solution  variables  (u,  11,0,  e)  at 
cell  centers  (Sullivan  et  ah,  2000,  2008)  which  leads  to  a  simple  compact  differencing  stencil. 
To  maintain  tight  velocity-pressure  coupling  the  flux  velocities  ( U,V )  are  located  at  cell  centers 
while  W  is  located  at  a  cell  face.  We  use  so-called  momentum  interpolation  to  construct  the 
intermediate  right  hand  sides  for  the  tendencies  of  the  flux  velocities.  The  spatial  differencing 
is  pseudospectral  in  the  horizontal  computational  directions  (£,  rj)  and  second-order  finite  differ¬ 
ences  in  the  ((-coordinate.  The  time  stepping  is  a  low-storage  third-order  Runge-Kutta  scheme, 
and  the  time  step  5t  is  picked  dynamically  based  on  a  fixed  Courant-Fredrichs-Lewy  (CFL) 
number. 

Compared  to  our  flat  LES  code,  the  major  algorithmic  change  is  the  pressure  formulation.  As 
a  consequence  of  the  incompressible  flow  assumption  and  the  non-orthogonal  mesh  we  are  forced 
to  solve  a  general  Poisson  equation.  A  direct  solution  of  this  elliptic  problem  is  not  possible,  and 
thus  we  use  the  stationary  iteration  scheme 


£>(1T+1)  =  £>(LP) 


1A 

5t  d 


Ufl  IT) 


(2) 


for  IT+1,  the  superscript  i  is  an  iteration  index  and  the  operator  T>  is  an  approximate  diago- 
nalization  of  the  full  Poisson  operator.  At  convergence  21(11*)  =  D(LL+1)  and  (2)  numerically 
satisfies  mass  conservation;  typically  20  to  30  iterations  are  required  to  drive  the  residual  error 
in  (2)  to  machine  accuracy  for  hard  problems.  (2)  is  designed  so  that  it  can  be  easily  inverted 
using  a  combination  of  2D  Fast  Fourier  Transforms  (FFTs)  and  tridiagonal  matrix  inversions, 
which  nicely  maps  onto  our  2D  Message  Passing  Interface  (MPI)  parallelization,  see  Fig.  5. 
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The  time  dependence  of  the  grid  modifies  the  LES  equations:  the  Jacobian  J(x,t)  appears 
inside  the  time  tendency  of  each  transport  equation  and  advection  contains  a  contribution  from 
the  grid  movement,  i.e.,  the  total  vertical  flux  of  variable  ip  depends  on  the  difference  {W  —  Zt)ip 
where  zt  is  the  grid  speed.  As  a  consequence  of  writing  the  equations  in  strong  conservation 
form,  we  determine  the  grid  speeds  based  on  the  discrete  solution  of  dJ/dt  =  dzt/d(,  which 
is  a  simplified  form  of  the  geometric  conservation  law  first  discussed  by  Thomas  and  Lombard 
(1979).  This  prevents  artificial  sources  and  sinks  from  developing  in  the  computational  domain. 

The  MPI  algorithm  utilizes  the  2D  brick  decomposition  shown  in  Fig.  5.  This  scheme  allows 
freedom  in  choosing  the  horizontal  and  vertical  decompositions  and  allows  simulations  with  a 
high  degree  of  spatial  anisotropy,  i.e.,  when  the  computational  box  is  wide  with  small  vertical 
extent.  The  algorithm  relies  on  matrix  transposes  to  carry  out  the  global  FFT  and  tri diagonal 
operations.  Further  details  about  the  LES  model  and  the  algorithm  are  given  by  Sullivan  et  al. 
(2010a, b);  Sullivan  and  Patton  (2011). 

2.  Numerical  experiments 

For  our  project  we  propose  carrying  out  the  7  simulations  listed  in  Table  1;  included  in  this 
table  is  an  estimate  of  the  computational  resources  required.  These  simulations  allow  us  to 
quantify  the  effects  of  thermal  stratification  and  a  change  in  wind  speed  on  the  turbulent  flow 
in  and  around  the  wave  field.  We  also  propose  a  stationary  bump  case,  i.e.,  the  wave  fields  will 
be  frozen  in  time  so  that  they  correspond  to  stationary  resolved  roughness  which  will  further 
expose  the  influence  of  wave  motion.  All  the  numerical  experiments  will  be  conducted  in  domains 
(Lx,  Ly,  Lz)  =  (1200,1200,800)  m  using  meshes  of  (Nx,  Ny,  Nz)  =  (1024,1024,512)  gridpoints. 
The  grid  spacing  in  the  horizontal  directions  is  constant  Ax  =  Ay  =  1.17  m  while  the  vertical 
grid  is  stretched  to  allow  fine  resolution  near  the  wavy  surface;  =  0.5  nr  or  less  at  the  surface. 
More  than  370  gridpoints  are  then  distributed  between  the  water  surface  and  the  boundary-layer 
inversion.  For  the  first  5  simulations  in  Table  1,  the  large  scale  winds  will  be  held  fixed  at  Ug  =  20 
ms-1  and  the  surface  buoyancy  will  be  varied  to  generate  cases  with  neutral,  unstable  and  stable 
stratification  to  match  the  conditions  in  the  Hi-Res  field  campaign  (Fig.  1).  Simulation  A10  will 
use  the  same  wavefields  as  N 20  but  the  large-scale  winds  will  be  reduced  to  10  ms-1  to  generate 
a  swell  dominated  condition,  a  commonly  observed  ocean  state.  Based  on  previous  simulations, 
the  impact  of  wave  age  and  swell  on  the  atmosphere  are  pronounced  as  shown  in  Fig.  4. 

The  computational  strategy  we  will  adopt  is  to  first  carry  out  simulations  with  a  flat  stationary 
lower  boundary  (F20,  FT0)  so  as  to  develop  near  equilibrium  3D  turbulence.  3D  volumes  from 
this  seed  simulation  will  then  be  used  to  spawn  simulations  with  different  combinations  of  wind, 
stratification,  and  waves.  This  strategy  makes  effective  use  of  the  computational  resources  since 
the  LES  code  with  a  flat  lower  boundary  runs  between  2  ~  3  times  faster  than  the  code  with 
moving  waves.  Hence  it  is  most  efficient  to  spinup  boundary-layer  turbulence  in  the  presence  of 
a  flat  lower  boundary. 

3.  Computational  resources  required 

To  estimate  the  core-hours  C  required  by  a  simulation  we  use  the  formula  C  =  Ncores  x  At  x 
A(stepS/3600.0  where  Ncore  is  the  number  of  cores,  At  is  the  core  seconds  required  per  compu¬ 
tational  step  and  Nstepi ,  is  the  number  of  steps  required  for  producing  equilibrium  turbulence. 
On  the  Cray  XE6  GarnetP  with  a  mesh  of  5122  x  128  using  512  cores  the  code  At  &  5.4  core- 
seconds/step.  The  number  of  time  steps  needed  to  achieve  6  large  eddy  turnover  times  is  105,  000 
and  thus  C  =  80,  640  core-hours.  To  estimate  the  computational  cost  for  the  much  higher  res¬ 
olution  simulations  on  Yellowstone  we  increase  the  Garnet  estimate  by  a  factor  of  39.5.  This 

^Details  about  Garnet  can  be  found  at  http://www.erdc.hpc.mil/hardware/index.htinl  . 
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Table  1:  Summary  of  proposed  simulations 


Run 

Comment 

Core-hour 

106 

Disk  space 
TB 

Archive  space 
TB 

N20 

neutral  stratification,  20  m/s 

1.6 

1.5 

4.5 

£20 

stable  stratification,  20  m/s 

1.6 

1.5 

4.5 

U20 

unstable  stratification,  20  m/s 

1.6 

1.5 

4.5 

B20 

neutral  stratification,  bumps 

1.6 

1.5 

4.5 

F20 

neutral  stratification,  flat  wall, 

20  m/s,  used  for  initial  conditions 

1.6 

1.5 

4.5 

N10 

neutral  stratification,  10  m/s 

1.6 

1.5 

4.5 

F 10 

neutral  stratification,  flat  wall, 

10  m/s,  used  for  initial  conditions 

1.6 

1.5 

4.5 

accounts  for  the  increase  in  resolution  (factor  of  19.75  that  includes  the  FFT  scaling  N  log  N) 
and  reduction  in  allowable  time  step  due  to  the  CFL  constraint  (factor  of  2).  We  then  reduce 
this  estimate  by  a  factor  of  2  assuming  Yellowstone  is  twice  as  fast  as  Garnet.  Thus  a  single 
simulation  on  Yellowstone  consumes  C  =  80,  640  x  39.5/2  =  1.6  x  106  core-hours.  The  total 
computational  cost  of  the  7  simulations  is  then  ~  12.0  x  106  core-hours. 

D.  Code  requirements  and  production-readiness 

1.  Code  to  be  used 

The  in-house  NCAR  wavy  LES  code  is  tested  and  production  ready  on  large  parallel  supercom¬ 
puters.  We  routinely  use  the  wavy  LES  code  on  the  three  major  supercomputer  platforms,  IBM, 
SGI,  and  Cray  all  with  different  compilers:  e.g.,  NCAR  IBM  ( Blueftre )  with  the  mpxlf90  com¬ 
piler;  ERDC  SGI  Altix  Ice  {Diamond)  with  the  ifort  compiler;  the  ERDC  Cray  XE6  ( Garnet ) 
with  the  pgf90  compiler.  The  LES  code  is  robust  across  different  platforms,  and  its  compliance 
to  the  F90  standard  and  the  use  of  the  most  common  functions  in  the  MPI  libraries  allows  the 
code  to  be  ported  to  new  machine  architectures  quickly.  A  flat  version  of  the  LES  code  (no  wavy 
boundaries)  was  one  of  the  benchmark  codes  used  in  the  procurement  for  the  NCAR  Wyoming 
Supercomputer  Center,  and  hence  we  expect  the  wavy  LES  code  will  port  to  Yellowstone  in  a 
straightforward  manner. 

2.  Programming  environment 

The  specifics  of  the  programming  environment  that  the  LES  code  utilizes  are  provided  in  Table 

2.  The  estimates  of  the  memory  usage,  wall-clock  time,  and  scratch  disk  space  are  based  on 
simulations  with  10242  x  512  «  540  x  106  gridpoints. 

3.  Benchmarks  and  scalability 

The  scaling  performance  of  the  flat  LES  code  for  three  different  machines  for  varying  workload 
as  a  function  of  the  total  number  of  processors  NP  is  documented  in  Sullivan  and  Patton  (2011). 
An  example  of  the  code  scaling  is  provided  in  Fig.  6  for  the  Cray  machines  at  the  National 
Energy  Research  Scientific  Computing  Center  (NERSC).  These  timing  tests  illustrate  the  present 
algorithm  exhibits  strong  scaling  over  a  wide  range  of  problem  sizes.  Note  the  code  is  able  to 
use  as  many  as  65,536  processors  for  a  large  problem  with  30723  gridpoints.  Generally,  the 
performance  only  begins  to  significantly  degrade  when  the  number  of  processors  exceeds  about  8 
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Table  2:  Programming  environment 


Programming  language: 
Parallel  paradigm: 

Estimated  cores/run: 
Memory  usage/core: 


F90  compliant 
MPI,  2D  domain 
decomposition 
4096  or  8192 
<  2  MB 


Required  libraries: 
Wall-clock  time/run: 
Scratch  disk  per  run: 
Restart  capability: 
Data  I/O: 


FFTPACK  and  MPI 
195  hours  on  8192  cores 
1.5  TB 

Yes,  routinely  used 
MPI  I/O 


times  the  smallest  grid  dimension  owing  to  increases  in  communication  overhead.  We  speculate 
that  the  performance  degradation  at  65,536  processors  is  a  consequence  of  the  small  node  memory 
on  Hopper  for  the  large  problem  tested,  30723  gridpoints.  The  code  also  exhibits  weak  scaling  ( 
i.e.,  where  the  problem  size  grows  as  the  number  of  processors  increases  so  the  amount  of  work 
per  processor  is  held  constant). 

Preliminary  tests  indicate  that  the  wavy  LES  code  scales  in  a  similar  manner  as  the  flat  code. 
The  wavy  LES  code  is  routinely  run  in  production  mode  with  either  512  or  1024  processors  on 
the  ERDC  machines  Diamond  and  Garnet  for  problems  with  33  x  106  gridpoints.  On  Yellowstone 
we  expect  to  routinely  use  4096  or  8192  processors  for  the  proposed  simulations  with  540  x  106 
gridpoints. 

E.  Data  analysis  and  visualization 

We  typically  use  Tecplot  360  (http://www.tecplot.com/)  on  our  desktop  machines.  This  visu¬ 
alization  package  allows  for  building  presentation  quality  animations  quickly  using  2D  planes  of 
data  with  time  dependent  moving  grids.  In  the  present  project  we  will  use  Tecplot’s  capabilities 
but  will  also  modify  our  analysis  to  use  NCAR’s  VAPOR  package  for  large  3D  visualizations. 
Also,  due  to  the  large  size  of  the  data  volumes  we  will  alter  our  usual  analysis  strategy,  leaving 
more  data  on  the  centralized  hie  and  data  storage  resources,  and  perform  extensive  analysis  and 
visualization  on  the  Geyser  and  Caldera  machines. 

F.  Data  management  plan 

We  estimate  that  our  simulations  require  approximately  5  ~  8  TB  of  working  disk  space  for  data 
analysis  and  visualization.  This  includes  2D  planes  sampled  at  an  adequate  temporal  resolution 
for  rapid  movie  making  of  selected  how  variables.  In  addition,  disk  space  is  needed  to  hold 
several  3D  volumes  for  conditional  sampling  and  non-routine  analysis  not  conducted  on  the  hy 
during  a  simulation.  Our  long  term  storage  plans  rely  on  keeping  several  archival  3D  volumes 
that  can  be  used  to  perform  code  restarts  if  further  analysis  is  needed.  The  LES  code  builds 
history  hies  with  numerous  statistics  as  function  of  time,  e.g.,  vertical  prohles  of  winds,  huxes, 
and  variances  in  addition  to  point  values.  These  history  hies  are  compact  and  provide  a  quick 
look  at  the  simulation  data  using  already  built  analysis  codes. 

G.  Accomplishment  report  on  prior  large  allocations 

The  Pis  have  worked  extensively  on  IBM,  SGI,  and  Cray  supercomputers  during  the  past  10 
years.  We  estimate  our  current  average  use  at  300,000  Gaus/year  on  NCAR’s  Bluefire,  more 
than  2  million  core-hours/year  at  NERSC  and  1  to  2  million  core-hours/year  on  DOD  ma¬ 
chines.  These  past  computational  experiments  focused  on  a  wide  range  of  LES  related  topics, 
viz.,  air-sea  interaction  in  the  marine  atmospheric  and  oceanic  boundary  layers,  turbulent  how  in 
vegetative  canopies  including  chemistry,  turbulent  dispersion,  how  over  small-scale  topography, 
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stable  boundary  layers  with  surface  heterogeneity,  and  effects  of  grid  resolution  on  LES  solutions. 
This  computational  work  has  produced  numerous  publications  that  appear  in  Journal  of  Fluid 
Mechanics ,  Journal  of  the  Atmospheric  Sciences,  Boundary-Layer  Meteorology,  Journal  of  Phys¬ 
ical  Oceanography,  Journal  of  Geophysical  Research,  and  Annual  Review  of  Fluid  Mechanics.  An 
itemized  list  of  recent  publications  can  be  found  at  http://www.mmm.ucar.edu/people/sullivan/ 
and  http://www.mmm.ucar.edu/people/patton/.  The  Pis  routinely  present  results  from  their 
computational  work  at  scientific  meetings  and  workshops.  The  Pis  research  is  funded  internally 
by  NSF  and  externally  by  the  Office  of  Naval  Research,  the  Army  Research  Office,  and  the 
Department  of  Energy.  The  Pis  regularly  collaborate  with  several  university  partners;  PPS  cur¬ 
rently  serves  on  5  PhD  thesis  committees,  viz.,  University  of  Colorado  (2),  Stanford  University 
(1),  UCLA  (1),  and  University  of  Rhode  Island  (1). 
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Figure  1:  Photograph  of  the  Floating  Instrument  Platform  (FLIP)  in  seas  generated  by  surface  winds 
of  15  to  20  ms _1  during  the  High  Resolution  Air-Sea  Interaction  (Hi-Res)  field  campaign  (2009)  off 
the  California  coast,  courtesy  E.  Terrill  Scripps  Institution  of  Oceanography.  The  goals  and  description 
of  Hi-Res  can  be  found  at  http://airsea.ucsd.edu/hires/.  Hi-Res  is  supported  by  the  Office  of  Naval 
Research. 
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Figure  2:  A  snapshot  of  the  wavefield 
height  h(x,  y,  t)  that  is  imposed  at  the 
bottom  of  the  LES  code  in  low  resolu¬ 
tion  runs,  h  is  built  from  a  sum  of  lin¬ 
ear  plane  waves.  Waves  propagate  left  to 
right  according  to  the  dispersion  relation¬ 
ship.  The  color  bar  is  in  units  of  meters. 


Figure  3:  An  instantaneous  x  —  z  slice  of  the  3D  time  varying 
computational  mesh  in  the  lowest  portion  of  the  atmospheric 
boundary  layer.  The  (£,  rj)  gridlines  become  level  surfaces  at 
about  100  m  above  the  water.  Only  a  fraction  of  the  grid  is 
displayed. 


Figure  4:  Snapshot  of  resolved  vertical  velocity  fluctuations  w' /u*  in  a  wave  following  x  —  y  plane 
near  the  water  surface  (  =  2.5  m.  The  left  panel  is  a  swell  dominated  regime  while  the  right  panel  is  a 
case  near  wind-wave  equilibrium.  The  wave  spectrum  at  the  bottom  of  the  PBLs  is  the  same.  Notice 
the  range  of  the  color  bar  is  different  between  the  two  cases. 


Figure  5:  An  example  of  2D  domain  decomposition  on  NP  =  9  processors.  The  number  of  horizontal 
processors  NPxy  =  3  and  the  vertical  domain  is  divided  into  NPZ  =  3  levels,  i.e. ,  the  total  number  of 
processors  NP  =  NPxy  ■  NPZ:  (a)  base  state  with  y  —  z  decomposition;  (b)  x  —  z  decomposition  used  for 
computation  of  y  derivatives  and  2D  planar  FFT;  and  (c)  x  —  y  decomposition  used  in  the  tridiagonal 
matrix  inversion  of  the  pressure  Poisson  equation. 
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Figure  6:  Performance  of  the  NCAR  parallel  LES  code  on  the  NERSC  Cray  supercomputers,  Franklin 
and  Hopper.  NP  is  the  total  number  of  processors  and  the  vertical  axis  is  total  computational  time 
t  X  NP  divided  by  total  work.  Nz  is  the  number  of  vertical  levels  and  MX}V  is  proportional  to  the  FFT 
work,  i.e.,  Mx  y  =  Nx>ylogNX!y  with  NXjV  the  number  of  gridpoints  in  the  x  and  y  directions.  Ideal  scaling 
corresponds  to  a  flat  line  with  increasing  number  of  processors.  This  figure  illustrates  that  the  LES  code 
exhibits  strong  scaling  for  different  combinations  of  problem  size  and  2D  domain  decompositions,  a) 
Green  ♦  problem  size  5123;  b)  Red  ♦  10243;  c)  Black  ♦  20483;  d)  Orange  ♦  30723;  and  e)  Blue  ♦  3072s. 
Cases  [a)  -  d)]  are  timing  tests  performed  on  Franklin  while  case  e)  is  a  timing  test  on  Hopper.  For  a 
given  number  of  total  processors  NP  the  symbols  are  different  vertical  and  horizontal  decompositions, 
i.e.,  different  combinations  of  (NPZ,  NPxy)  as  depicted  in  Fig.  5. 
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